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ABSTRACT 

We present a new semi-analytic model that self-consistently traces the growth of 
supermassive black holes (BH) and their host galaxies within the context of the 
ACDM cosmological framework. In our model, the energy emitted by accreting 
black holes regulates the growth of the black holes themselves, drives galactic 
scale winds that can remove cold gas from galaxies, and produces powerful jets 
that heat the hot gas atmospheres surrounding groups and clusters. We present 
a comprehensive comparison of our model predictions with observational mea- 
surements of key physical properties of low-redshift galaxies, such as cold gas 
fractions, stellar metallicities and ages, and specific star formation rates. We 
find that our new models successfully reproduce the exponential cutoff in the 
stellar mass function and the stellar and cold gas mass densities at z 0, and 
predict that star formation should be largely, but not entirely, quenched in mas- 
sive galaxies at the present day. We also find that our model of self-regulated BH 
growth naturally reproduces the observed relation between BH mass and bulge 
mass. We explore the global formation history of galaxies in our models, present- 
ing predictions for the cosmic histories of star formation, stellar mass assembly, 
cold gas, and metals. We find that models assuming the "concordance" ACDM 
cosmology overproduce star formation and stellar mass at high redshift (z ^ 2). 
A model with less small-scale power predicts less star formation at high redshift, 
and excellent agreement with the observed stellar mass assembly history, but 
may have difficulty accounting for the cold gas in quasar absorption systems at 
high redshift {z ~ 3-4). 
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1 INTRODUCTION 

It is now well-established that the Cold Dark Mat- 
ter p aradigm for structure formation l|Blumenthal et al] 
Il984l ). in its modern, dark energy-dominated (ACDM) 
incarnation, provides a remarkably successful paradigm 
for interpreting a wide variety of observations, from the 
cosmic micr owave background fluctuations at z 1000 
l|Spergel et al.. .2003. . .2007. ') , to the large-scale clustering 



of galaxies at z ~ (iPercival et al.ll2002l : iTegmark et al.l 
[2004; Eiscn stein et al.ll2005l ). Other successful predictions 
of the ACDM paradigm include the cos mic shear field as 
measured by weak gravitational lensing (iHeymans et al.l 
I2OO5I : iBacon et al]|2005l : iHoekstra et al.ll2006l ). the small- 



scale power spectru m as probed b y the Lyma n- 



forest (Dcsiacqucs fc Nusserl l2005l: iJena et abl 



and the number densities of l|Borgani et al.l | 200lf). and 



alph a 



20051 ). 



baryon fractions withi n, galaxy clusters l|Allen et al.l 
I2OO4I : IWhite et al.lll993D . 



E-mail: somerviIIe@mpia.de 



ACDM as a paradigm for understanding and sim- 



© 0000 RAS 



2 Somerville et al. 



ulating galaxy formation has had more mixed success. 
In thi s picture, as originally proposed by IWhite fc Re3 
l| 19781 ) and [Blumenthal et al. (1984), galaxies form when 
gas cools and condenses at the centers of dark-matter 
dominated potential wells, or "halos". More detailed 
calculations, using semi-analytic and numerical sim- 
ulations of galaxy formation, have shown that this 
framework does provide a promising qualitative under- 
standing of many features of galax ie s and their evo- 
lution (e.g.. iKauffmann et al.l il993l: ICole et al.l 1 1994 : 



Somerville fc Primackl Il999l: iKauffmann et al.| 11998 
Cole et al.1 I2OO0I : ISomerville et al.1 I2001I : iBaughl I2OO6I ') 



However, it has been clear for at least a decade now that 
there is a fundamental tension between certain basic pre- 
dictions of the ACDM-based galaxy formation paradigm 
and some of the most fundamental observable proper- 
ties of galaxies. In this paper we focus on two intercon- 
nected, but possibly distinct, problems: 1) the "overcool- 
ing" or "massive galaxy" problem and 2) the star forma- 
tion "quenching" problem. 

The first problem is manifested by the fact that both 
semi-analytic and numerical simulations predict that a 
large fraction of the available baryons in the Universe 
rapidly cools and condenses, in conflict with observations 
which indicate that only about ~ 5-10 % of the baryons 
are in the form of cold g as and stars l|Bell et al ] |2003al : 
iFukugita fc Peeblei2004l ). This is due to the fact that gas 
at the densities and temperatures characteristic of dark 
matter halos is expected to cool rapid ly, and is related 
to the classical "c ooling flow" problem (Fabian & Nulsen' 
ll977l : ICowie fc Binnev 1977: Mathews & Brcgrnan 1978) . 
Direct observations of the X-ray properties of hot gas in 
clusters similarly imply that this gas should have short 
cooling times, particularly near the center of the clus- 
ter, but the condensations of stars and cold gas at the 
centers of these clusters are much smaller than would be 
expected if the hot gas had been cooling so e fflciently over 
the lifetime of the cluste r (for reviews see lFabianlll993 : 
iPeterson fc Fabian|[200^ ). Moreover, X-ray spectroscopy 
shows that very little gas is cooling below a temperature 
of about one-third of the virial temperature of the cluster 
l|Peterson et al.ll2003D . 

A further difflculty is that there is a fundamen- 
tal mismatch between the shape of the dark matter 
halo mass function and that of the observed mass func- 
tion of cold baryons ( cold g as and stars) in g alaxies 
(|SomerviUe fc Primackl Il999l : iBenson et all I2OO3I I. The 
galaxy mass function has a sharp exponential cutoff 
above a mass of about a few times 10^'^Mq , while the halo 
mass function has a shallower, power-law cutoff at much 
higher mass (~ few x lO^^M©). There is also a mismatch 
at the small-mass end, as the mass function of dark mat- 
ter halos is much steeper than that of galaxies. If we are 
to assume that each dark matter halo hosts a galaxy, this 
then implies that the ratio between the luminosity or stel- 
lar mass of a galaxy and the mass of its dark matter halo 
varies strongly and non-monotonical l y with halo mass 
jKra-vtsov et al.1 l 2004a: Conro v et al] I2OO6I : |Wang ot ajj 
I2OO6I : iMoster et al.. ,2008). such that "galaxy formation" 
is much more inefflcient in both small mass and large 
mass halos, with a peak in efficiency close to the mass of 
our own Galaxy, ~ lO^'^M©. On very small mass scales 



(below halo velocities of ~ 30 — 50 km s~ ), the col- 
lapse and cooling of baryons may be suppressed by the 
presence of a photoionizing background (iEfst athjou 199^; 
iThoul fc Weinberg! Il996l : iQuinn et al]|l996l ). For larger 
mass halos (up to Vvir — 150 — 200 km s~^), the standard 
assumption is that winds driven by massive stars and su- 
pernovae are able to heat and expell gas, resulting in low 
baryon fractions in small mass halos (.Wh ite & Roes 19781: 
iDekel fc Sildll986l : IWhite fc FrenkHiggii ). However, stel- 
lar feedback probably cannot provide a via ble solution to 
the o vercooling problem in massive halos (jBenson et al.l 
I2OO3I ): stars do not produce enough energy to expell 
gas from these large potential wells; and the massive, 
early type galaxies in which the energy source is needed 
have predominantly old stellar populations and little or 
no ongoing or recent star formation. Other solutions, 
like thermal conduction, have been explored, but prob- 
ably do not provide a full solution (|Benson et al.ll2003l : 
IVoigt fc Fab"ianll2004l ). 

The second problem is related to the correlation 
of galaxy structural properties (morphology) and spec- 
trophotometric properties (stellar populations) with stel- 
lar mass. The presence of such a correlation has long been 
known, in the sense that more massive galaxies tend to be 
predominantly spheroid-dominated, with red colors, old 
stellar populations, low gas fractions, and little recent 
star formation, while low-mass galaxies tend to be disk- 
dominated and gas-rich, with blue colors and ongoing star 
formation ('e.g. [Roberts fc Havnes|[l994l ). More recently, 
with the advent of large galaxy surveys such as the Sloan 
Digital Sky Survey (SDSS), we have learned that the 
galaxy color distribution (and t hat of other related prop- 
erties) is strongly bimodal (e.g. iBaldrv et al]|2004l '). and 
that the transition in galaxy properties from star-forming 
disks to "dead" spheroids occurs rather sharp ly, at a char- 
acteristic stellar mass of ^ 3x 1 0"Mq (jKauffmann et al.l 
l2003al : fBrinchmann et al.ll2004l ). In contrast, the "stan- 
dard" ACDM-based galaxy formation models predict 
that massive halos have been assembled relatively re- 
cently, and should contain an ample supply of new fuel for 
star formation. These models predict an inverted color- 
mass and morphology-mass relation (massive galaxies 
tend to be blue and disk dominated, rather than red and 
spheroid dominated) and no sharp transition or strong 
bimodality. Thus, the standard paradigm of galaxy for- 
mation does not provide a physical explanation for the 
"special" mass scale (a halo mass of ~ IO^'^Mq, or a stel- 
lar mass of 3 X 10^" Mq) which marks both the peak of 
galaxy formation efficiency and the transition in galaxy 
properties seen in observations. 

Several pieces of observational evidence provide clues 
to the solution to these problems. It is now widely be- 
lieved that every spheroid-dominated galaxy hosts a nu- 
clear supermassive black hole (SMBH), and that the 
mass of the SMBH is tightly correlated with the lu- 
minosity, mass, and velocity dis persion of the stellar 
spheroid (iKormendv fc Richston^ 1995; Ma gorrian et al.l 
19981: iFerrar ese fc Merritt' '2OOO': Gcbh ardt et al.1 I2OO0I : 



Marconi fc Hunt ,2003 : Haring fc Rix,2004r ). These cor- 
relations may be seen as a kind of "fossil" evidence that 
black holes were responsible for regulating the growth of 
galaxies or vice versa. This also implies that the most 
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massive galaxies, where quenching is observed to be the 
most efficient, host the largest black holes, and there- 
fore the available energy budget is greatest in precisely 
the systems where it is needed, in contrast to the case 
of stellar feedback. The integrated energy released over 
the lifetime of a SMBH (~ 10^" - 10*^^ erg) is clearly 
very sig nificant compared with galaxy binding energies 
(|SUk fc R ocs 1998). In view of these facts, it seems al- 
most inconceivable that AGN feedback is not important 
in shaping galaxy properties. 

However, in order to build a complete, self-consistent 
machinery to describe the formation and growth of black 
holes within the framework of a cosmological galaxy for- 
mation model, and to attempt to treat the impact of the 
energy feedback from black holes in this context, we need 
to address several basic questions: 1) When, where, and 
with what masses do seed black holes form? 2) What trig- 
gers black hole accretion, what determines the efficiency 
of this accretion, and what shuts it off? 3) In what form 
is the energy produced by the black hole released, and 
how does this energy couple with the host galaxy and 
its surroundings? In order to address some of these ques- 
tions, we first identify two modes of AGN activity which 
have different observational manifestations, probably cor- 
respond to different accretion mechanisms, and have dif- 
ferent physical channels of interaction with galaxies. 



1.1 The Bright Mode of Black Hole Growth 

Classical luminous quasars and their less powerful 
cousins, optical or X-ray bright AGN, radiate at a signifi- 
cant fraction of th e ir Eddington limit (L ~ (0.1 — l)Z/Edd; 
IVestergaardI l2004l : iKoUmeier |2006| ). and are be- 

lieved to be fe d by optically thick, g eometrically thin ac- 
cretion disks l|Shakura fc Svunvaevlil973 '). We will refer 
to this mode of accretion as the "bright mode" because 
of its relatively high radiative efficiency (with a fraction 
f?rad ~ 0.1 — 0.3 of the accreted mass converted to radi- 
ation). The observed space density of these quasars and 
AGN is low compared to that of galaxies, implying that 
if most galaxies indeed host a SMBH, this "bright mode" 
of accretion is only "on" a relatively small fraction of 
the time. Constraints from quasar clustering and vari- 
ability imply that quasar l i fetimes must be ^ 10^'^ y r 
(|Martini fc Weinbergj l200ll : iMartini fc Schneider] 120031 '). 
These short timescales combined with the large observed 
luminosities immediately imply that fueling these objects 
requires funneling a quantity of gas comparable to the en- 
tire supply of a large galaxy (~ 10^ — lO^^M©) into the 
central regions on a timescale of order the dynamical time 
(~ few X 10^ - 10* yr). 

These considerations alone lead one to consider 
galaxy-galaxy mergers as a promising mechanism for 
triggering this efficient accretion onto nuclear black 
holes. The observational association of mergers with 
enhanced star formation, particularly with the most 
violent observed episodes of star formation exhibited 
by Ultra Lumi nous Infrared Galaxies ( ULIRGS), is 
well-established llSanders fc Mirabell [l996l: iFarrah et al, 
200l':'Colina et al."200l': 'Barton et al."2000': 'Woods et al. 



l2007l : iLi et all [2007l 'l. Moreover, numerical simulations 
have shown that tidal torques during galaxy mergers 
can drive the rapid inflows of gas that are needed 
to fuel both the intense starbursts and rapid black 
hole accretion associated w i th ULIRGS and quasars 
' Hernauistlll989l:lBarneslll992l:lBarnes fc Hernauisdll996| : 



Mihos fc He rnguist '199 
Pi Matteo et a l. 20oi) 



19961 : ISpringel et al.l l2005bl : 

As well, it seems that if one 
can probe sufficiently deep to study the SED be- 
neath the glare of the quasar, one always uncovers ev- 
idence of y oung stellar popu l ation s indicative of a recent 
starb ur st (iBrotherton et al.l 199"p : [Canalizo fc StocktorJ 

il 



20011: 



Sanchez et al. 



Kauffmann et all l2003bl: iJahnke et al.. ,2004 
' 20041 : IVanden Berk et al.l 120061 ). Near- 
equal mass (major) mergers also have the attractive fea- 
ture that they scramble stars from circular to random or- 
bits, leading to morphological transformation from disk 
to spheroid (iToomre fc Toomre|[l973 : lBarneslll988l , Il992l : 
lHernauistlll992l . 19931 ). If spheroids and black holes both 
arise from violent mergers, this provides a possible ex- 
planation for why black hole properties always seem to 
be closely associated with the spheroidal components of 
galaxies. 

What impact does the energy associated with this 
rapid, bright mode growth have on the galaxy and on 
the growth of the black hole itself? Long thought to 
be associated only with a small subset of objects (e.g. 
Broad Absorption Line (BAL) quasars), high- velocity 
winds have been detected in a variety of different types of 
quasar syste ms (dc Kool et al. 2001; Po unds et al.l 120031 : 
IChartas et al . 2003; Pounds & P age 200i), and are now 
believed to be quite ubiquitous jCangulv fc BrothertonI 
I2OO8I ). However, their impact on the host galaxy remains 
unclear, as the mass outflow rat es of these winds are dif- 
flcult to constrain ( though see Steenbrugge et al.l [20051 : 
IChartas et aL 20"07l : iKrongold et al.ll2007l ). Recently, nu- 
merical simulations of galaxy mergers including black 
hole growth found that depositing even a small fraction 
(~ 5 %) of the energy radiated by the BE into the ISM 
can not only halt t he accretion onto the BH , but can drive 
large-scale winds l|Di Matteo et all |2005| ) . These winds 
sweep the galaxy nearly clean of cold gas and halt fur- 
ther star formation, leaving behind a rapi dly reddening, 
spheroidal remnant (Springol et al. .2005al ). 

To study how the interplay between feedback from 
supermassive black hole accretion and supernovae, galaxy 
structure, orbital conflguration, and gas dissipation com- 
bine to determine the properties of spheroidal galax- 
ies formed through mergers, hundreds of hydrodynam- 
i cal simulations wer e performed b y Robertson et al. 
(|2006bl : 120063 : l2006al ) and Cox e t al. I2006bl;l2006al )^ 



2006; Woods fc Gcllcr 2007; Barton et al. .2007.;.Lin et al. 



the methodology p r esented by iDi Matteo et al.) l|200oi ) 
and ISpringel et al] (l2005d ). Robertson et al. and Cox 
et al. analyzed the merger remnants to study the red- 
shift evolution of the BH mass-a relation, the Funda- 
mental Plane, phase-space density, and kinematic prop- 
erties. This extensive suite has been supple mented by 
additional simulations of minor mergers from I Cox et al] 
(|200S ). Throughout the rest of this paper, when we refer 
to "the merger simulations", we refer to this suite. 

Based on their analysis of these simulations, 
iHopkins et al.l (|2007al ) have outlined an evolutionary se- 
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quence from galaxy-galaxy merger, to dust-enshrouded 
starburst and buried AGN, blow-out of the dust and 
ISM by the quasar- and starburst-driven winds, to clas- 
sical (unobscured) quas ar, post-starbu r st gala xy, and fi- 
nally "dead" elliptical. iHopkins et~aLl (l2007bl ) find that 
in the merger simulations, the accretion onto the BH 
is eventually halted by a pressure-driven outfiow. Be- 
cause the depth of the spheroid's potential well deter- 
mines the amoun t of momentum neces sary to entrain 
the infalling gas, iHopkins et al.l l|2007bl ) find that this 
leads to a "Black Hole Fundamental Plane", a corre- 
lation between the final black hole mass and sets of 
spheroid structural/dynamical properties (mass, size, ve- 
locity dispersion) sim i lar to the one seen in o bservations 
l|Hopkins et"aLll2007d : [Marconi fc H unt"2003'). 

Furthermore, Hopkins et al. ((2005d; 2005i; l2005bl ) 
have shown that the self-regulated nature of black hole 
growth in these simulations leads to a characteristic form 
for quasar lightcurves. As the galaxies near their final 
coalescence, the accretion rises to approximately the Ed- 
dington rate. After the critical black hole mass is reached 
and the outflow phase begins, the accretion rate enters a 
power-law decline phase. Although most of their growth 
occurs in the near-Eddington phase, quasars spend much 
of their time in the decline phase, and this implies that 
many observed low-luminosity quasars are actually rela- 
tively m assive black ho l es in th e last stages of their slow 
decline. iHopkins et al] l|2006dh found that when these 
lightcurves are convolved with the observed mass func- 
tion of merging galaxies, the predicted AGN luminosity 
function is consistent with obse r vation s. Moreover, Hop- 
kins et al. l|2005bl : l2005d : l2006al : l2006d ) have shown that 
this picture reproduces many quasar and galaxy observ- 
ables that are difficult to account for with more simpli- 
fied assumptions about QSO lightcurves, such as differ- 
ences in the quasar luminosity function in different bands 
and redshifts, Eddington ratio and column density distri- 
butions, the X-ray background spectrum, and relic red, 
early type galaxy population colors and distributions. 



1.2 The Radio Mode 

The second mode of AGN activity is much more com- 
mon, and in general less dramatic. A fairly large frac- 
tion of massive galaxies (particularly galaxies near the 
centers of g roups and clusters) a re detected at radio 
wavelengths (|Best et al.ll2005l . l2007t ). Most of these radio 
sources do not have emission lines char acteristic of clas- 
sical optical or X-ray bright quasars (|Best et al.l l2005l : 
iKauffmann et all |2007| ). and their accretion rates are 
believed to be a sm all fraction of the Eddington rate 
l|Raffertv et al.' 2006*) . They are extremely radiatively in- 
efficient jBirz an ct al. 2004), and thought to be fuelled by 
optically thin, geometrically thick accretion as expected 
in ADAF and ADIO S m odels such as those prop o sed b y 
iNaravan fc yH (|1994| ) and lBlandford fc BegelmanI (|l999l ) . 
Because these objects are generally identified via their 
radio emission, we refer to this mode of accretion and 
BH g rowth as the "radio mode" f following ICroton et al.l 
12009 ). 

Although these black holes seem to be inefficient at 



producing radiation, they can apparently be quite ef- 
ficient at producing kinetic energy in the form of rel- 
ativistic jets. Intriguingly, the majority of cooling flow 
clusters host these active ra dio galaxies at their centers 
(iDunn fc Fabianll2006l . l2008l ). and X-ray maps reveal that 
the radio lobes are often spatially coincident with cavi- 
ties, thought to be bubbles flU ed with rel ativistic plasma 
and inflated by the jets (Mc Namara fc Nulsen 2007. . and 
references therein) . The observations of these bubbles can 
be used to estimate the work required to inflate them 
against the pressure of the hot medium ( Birzan et al.l 
l2004l : iRaffertv et alll2006l : lAUen et al. l l2oodr and hence 
obtain lower limits on the jet power. 

While the idea that radio jets provide a heat 
source that could counteract c ooling flows has been 
discussed for many y ears fe.g. iBinnev fc Taborl Il995l: 
Churazov et al.l |2002| : iFabian et al.l |2003| : lOmma et al.l 



2004l : lBinnevll2004h . these observations now make it pos- 



sible to investigate more quantitatively whether the heat- 
ing rates are sufficient to offset the cooling rates in 
groups and clusters. Several studies conclude that in 
the majority of the systems studied, the AGN heat- 
ing traced by the power in the X-ray cavities alone 
is comparable to or in exces s of the energy being 
ra diated by the cooling gas ^McNamara ct al.' 20061; 
i Raffertv et al.ll2006l: iFabian et al. 2006: B est ct al. 20oi; 
iMcNamara fc Nulsenll2007l : lDunn fc FabianlboOSt ). More- 
over, the net cooling rate is correlated with the observed 
star formation rate in the central cD galaxy, indicating 
that ther e may be a self-regulating cycle of heating and 
cooling (iRaffertv et al.ll200Q ). 

Several other physical processes that could sup- 
press cooling in large mass halos have been sug- 
:ested and explored, such as th ermal conduction 



■Benson et a l. 2003: 'Voigt fc Fabianll2004 ). multi-phase 
coohng ifMaller fc Bullock ,2003 ). or heating by sub- 
struc t ure or clumpy accretion (iKhochfar fc Ostrikerl 
l2007l : iNaab et al.ll2007l : iDekel fc Birnboimll2008l ). While 
some or all of these processes may well be important, in 
this paper we will investigate whether it is plausible that 
"radio mode" heating alone can do the job. 



1.3 A Unified Model for Black Hole Activity 
and AGN Feedback 

All of this begs the question: what determines whether a 
black hole accretes in the "bright mode" or "radio mode" 
state? An interesting possi ble answer comes from an 
analogy with X-ray binaries (|Jesterll2005l : iKording et al.l 
I2OO6I ) . Observers can watch X-ray binaries in real time as 
they transition between two states: the "low/hard" state, 
in which a steady radio jet is present and a hard X-ray 
spectrum is observed, and the "high/soft" state, in which 
the jet dissapears a nd the X-ray spectrum has a soft, 
therm al component l|Maccarone et al ] |2003l : lFende7 et al.l 
I2OO4I ). The transition between the two states is thought to 
be connected to the accretion rate itself: the "high/soft" 
state is associated with accretion rates of Si (0.01-0.02) 
TTiEdd and the existence of a classical thin accretion disk, 
while the "low/hard" state is associated with lower ac- 
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cretion rates and radi ative^ inefficient ADAF/ADIOS 
accretion ('Fender et al." 2004). 

Recently, Siiacki ct al. (200'^) have applied this idea 
in cosmological hydrodynamic simulations, by assuming 
that when the accretion rate exceeds a critical value, 
"bright mode" feedback occurs (AGN-driven winds), 
while when the accretion rate is lower, "radio mode" (me- 
chanical bubble feedback) is implemented. The results 
of their simulations appear promising — they produced 
black hole and stellar mass densities in broad agreement 
with observations. In addition, they found that their 
implementation of AGN feedback was able to suppress 
strong cooling flows and produce shallower entropy pro- 
files in clusters, and to quench star formation in massive 
galaxies. However, the very large dynamic range required 
to treat the growth of black holes and galaxies in a cos- 
mological context — from the sub-pc scales of the BIf ac- 
cretion disk to the super-Mpc scales of large scale struc- 
ture — means that numerical techniques such as these 
will likely need to be supplemented by semi-analytic or 
sub-grid methods for some time to come. 

Our approa ch is in many resp ects very similar in 
spirit to that of ISiiacki et al.l (|2007l ). although of course 
we are forced to implement both modes of AGN feedback 
in an even more schematic manner because we are using a 
semi-analytic model rather than a numerical simulation. 
We adopt fairly standard semi-analytic treatments of the 
growth of dark matter halos via accretion and mergers, 
radiative cooling of gas, star formation, supernova feed- 
back, and chemical evolution. We then adopt the picture 
of self-regulated black hole growth and bright mode feed- 
back in mergers discussed in ijl.ll and implement these 
processes in our model using the results extracted from 
the merger simulations described above. We assume that 
the radio mode is fueled instead by hot gas in quasi- 
hydrostatic halos, and that the accretion rate is described 
by Bondi a ccretion from an i sothe rmal cooling flow as 
proposed bv lNulsen fc FabianI (|2000l . NFOO). We cahbrate 
the heating efficiency of the associated radio jets against 
direct observations of bubble energetics in clusters. 

A number of authors have previously explored the 
formation of black holes and AGN in the context of 
CDM-basod semi-analytic models of varying complex- 
ity l|Efstathiou fc Ro cs 1988; Ka uffmann fc I fachnclt 
20001: IWvithe fc Loebl [20021: feromlev et al.1 [2004 



Scannapieco fc OhI |2004| : IVolonteri et al.1 l2003l : 
Volonteri fc Reej2005h . and recently several studies have 



also investigated the impact of A GN feedback on galaxy 
formation using such models (ICattaneo et al] 
Croton et aLll2006l : [Bower et al.ll2006l : iMenci et al 



200e 



200e 



Schawinski et al.[ [20061 : [Kang et al.[ [20061 : [ Monaco et al.l 
20071 ). The models that we present here differ from previ- 
ous studies of which we are aware, in two main respects: 
1) we implement detailed modelling of self-regulated 
black hole growth and bright mode feedback based on 
an extensive suite of numerical simulations of galaxy 
mergers and 2) we adopt a simple but physical model 
for radio mode accretion and heating, and calibrate 
our model against direct observations of accretion rates 
and radio jet heating efficiencies. We present a broader 
and more detailed comparison with observations than 
previous works, and highlight some remaining problems 



that have not previously been emphasized. As well, 
unlike most previous studies, we calibrate our models 
and make our comparisons in terms of "physical" galaxy 
properties such as stellar mass and star formation rate, 
which can be estimated from observations, rather than 
casting our results in terms of observable properties such 
as luminosities and colors. Our results are therefore less 
sensitive to the details of dust and stellar population 
modelling, and easier to interpret in physical terms. 

The goals of this paper are to present our new models 
in detail, and to test and document the extent to which 
they reproduce basic galaxy observations at z = and 
the global cosmic histories of the main baryonic compo- 
nents of the Universe. The structure of the rest of this 
paper is as follows. In ^ we describe the ingredients of 
our models and provide a table of all of the model pa- 
rameters. In ^ we present predictions for key properties 
of galaxies at z ~ and for the global history of the 
main baryonic components of the Universe: star forma- 
tion, evolved stars, cold gas, metals, and black holes. We 
conclude in §31 



2 MODEL 

Our model is based on the semi-analytic ga l axy fo rmation 
cod e described in [Soinervifie fc Primack[ ([l999l . SP99) 
and [Somerville et al.l l|200ll . SPFOl), with several major 
updates and important new ingredients, which we de- 
scribe in detail here. Unless specified otherwise, we adopt 
the "Concordance" ACDM model (C-ACDM), with the 
parameters given in Table [1] which has been used in 
many recent semi-analytic studies of galaxy formation. 
In i)3.5l we also consider a model that uses the set of 
parameters obtained from the three year results of the 
Wilkinson Microwave Anisotropy Probe bv lSpergel et al.l 
(|2007[ ). which are also specified in Table[T] We refer to this 
as the "WMAP3" model. We assum e a universal Ch abrier 
stellar initial mass function (IMF: [Chabrie3 [20031 ). and 
where necessary we convert all observations used in our 
comparisons to be consistent with this IMF. 



2.1 Dark Matter Halos, Merger Trees and 
Substructure 

We compute the number of "root" dark matter (DM) ha- 
los as a functio n of mass at a desired ou tput redshift us- 
ing the model o f [Sheth fc Tormenl(|l999[ ). which has been 
shown to agree well with numerical simulations. Then, for 
each "root" halo of a given mass Mo and at a given output 
redshift, we construct a realization of the merger history 
based on the method described in [Somerville fc Kolatt[ 
(|l999l . SK99). We have introduced a modification to the 
SK99 algorithm, which we find leads to better agreement 
with N-body simulations. We choose the timestep At by 
requiring that the average number of progenitors Np be 
close to two, by inverting the equation for Afprog(Mo, At) 
(see SK99). We then select progenitors as described in 
SK99, but do not allow the number of progenitors to ex- 
ceed Np + \J Np + 1. We follow halo merging histories 
down to a minimum progenitor mass of IO'^'^Mq, and 
our smallest "root" halos have a mass of 10" Mp,. We 
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Table 1. Summary of Cosmological Parameters 



parameter 



description 



Concordance ACDM WMAP3 



cosmological parameters 

VIA 

Ho 

fb 

Us 



present day matter density 0.30 

cosmological constant 0.70 

Hubble Parameter [km/s/Mpc] 70.0 

cosmic baryon fraction 0.14 

power spectrum normalization 0.9 

slope of primordial power spectrum 1.0 



0.2383 
0.7617 

73.2 
0.1746 
0.761 
0.958 



Table 2. Summary of the galaxy formation parameters in our "fiducial" model. We also specify the section in the paper where 
a more detailed definition of each set of parameters can be found, and whether the parameter is considered to be fixed based on 
direct observations or numerical simulations (F), or adjusted to match observations (A). 



parameter 


description 


fiducial value 


fixed / adj usted 


photoionization squelching (92.3D 

^overlap i ^reionize 


rcdshift of ovcrlap/rcionization 


11, 10 


F 


Quiescent star formation ff|2.5.1D 

^Kcnn 
AfK 
Xgaa 
^crit 


normalization of Kennicutt Law [MQyr~-^ kpc~^] 
power law index in Kennicutt Law 
scale radius of gas disk, relative to stellar disk 
critical surface density for star formation [M0pc~^] 


8.33 X 10-5 
1.4 
1.5 
6.0 


A 
F 

A 

A 


burst star formation fQ2.5.2D 

A'crit 
^burst 
Tburst 
'Haurst 


critical mass ratio for burst activity 
burst efficiency for 1:1 merger 
dependence of burst efficiency on mass ratio 
burst timescale 


0.1 
eqn. |9] 
eqn.E] 
eqn. [TO] 


F 
F 
F 
F 


merger remnants morphology (92. 6D 

/sph 
/scatter 


fraction of stars in spheroidal remnant 
fraction of scattered satellite stars 


ean.im 
0.4 


A 
A 


supernova feedback (^2.7\t 

^ject 
Xreinfall 


normalization of reheating function 
power law slope of reheating function 
velocity scale for ejection of reheated gas [km/s] 
timescale for re-infall of ejected gas 


1.3 
2.0 
120 
0.1 


A 
A 
A 
A 


chemical evolution (92.8D 

y 

R 


chemical yield (solar units) 
recycled fraction 


1.5 
0.43 


A 
F 


black hole growth (^2^9)1 

'?rad 

/bH, final 
/BH.crit 


efficiency of conversion of rest mass to radiation 
mass of seed BH [Mq] 

scaling factor for mass of BH at end of merger 
scaling factor for "critical mass" of BH 


0.1 
100 
2.0 
0.4 


F 
F 
A 
F 


AGN-driven winds Ca2.10ll 

^wind 


effective coupling factor for AGN driven winds 


0.5 


F 


radio mode feedback (92.1111 

'^radio 
^hcat 


normalization of "radio mode" BH accretion rate 
coupling efficiency of radio jets with hot gas 


3.5 X 10^3 
1.0 


A 

F 



have also implemented our models within N-body based 
merger trees, and do not find any significant changes to 
our results. 

We assign two basic properties to every dark mat- 
ter halo in each of our merger trees: the angular momen- 
tum or spin parameter, and the concentration parameter. 



which describes the matter density profile. We express the 
angular momentum in terms of the di mensionless spi n 
parameter A = AlBftl^/^G^^Mvir"''^^ (iPeebled [l969l ). 
where Eh is the total energy of the halo and Mvir is the 
virial mass. Numerical N-body simulations have demon- 
strated that A is uncorrelated with the halo's mass and 
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concentration (|Bullock et all l2001al : iMaccio etall [2003) 
and does not evolve with redshift. The distribution of A 
is log-normal, with m ean A = 0.05 and width a\ = 0.5 
l|Bullock et al.ll2001al ). We assign each top-level halo a 
value of A by selecting values randomly from this distri- 
bution, assuming that it is not correlated with any other 
halo properties or with redshift. The halo at the next 
stage of the merger tree inherits the spin parameter of its 
largest progenitor. 

We assume that the initial density profile of each 
halo i s described by the Navarro-Frenk- White (NFW) 
form ([Navarro et al.' 1997"), and compute the character- 
istic concentration parameter cnfw for the appropriate 
mass and redshift u sing a fitting formul a based on nu- 
merical simulations (iBuUock et al.ll2001bh . We adopt the 
updat ed normalization of cnfw (Mvir) from[M accio et al.l 
([2003). We neglect the scatter in cnfw at fixed mass, 
as well as the k nown correlation betw een cnfw and halo 
merger history (|Wechsler et al.ll2002l ). 

At each stage in the merging hierarchy, one or more 
halos merge together to form a new, virialized dark mat- 
ter halo. The merged halos (hereafter referred to as "sub- 
halos") and their galaxies, however, can survive and con- 
tinue to orbit within the potential well of the parent 
DM halo for some time. The time it takes for the satel- 
lite to lose all of its angular momentum due to dynam- 
ical friction and merge with the central galaxy is typi- 
cally modelled with some variant of the Chandrasekha r 
formula (see e.g. §2.8 of ISomerville fc Primackl Il999l ). 
Here, we use an updated v ersion of this formula from 
iBovl an-Kolchin et al] (|2008l ). which accounts for the tidal 
mass loss of sub-halos as they orbit within the host halo, 
as well as the dependence on the energy and angular mo- 
mentum of the orbit. Because the merger time is pro- 
portional to Mhost/A^sat, accounting for this mass loss 
increases the time it takes for small mass satellites to 
merge. 

These satellites may eventually lose so much of their 
mass tha t they become tidally disru pted. Based on the 
result s of lTavlor fc Babull \2QQ4 ) and lZentner fc BuUockl 
(|2003l ). we assume that satellites lose ~ 30 — 40 percent of 
their mass per orbital period, and that when the mass has 
been stripped down to the mass within the NFW scale 
radius = rvir/cNFW, we consider the satellite to be 
tidally destroyed. Sub-halos that survive until they reach 
the center of the parent halo are assumed to merge with 
the central object. Subhalos that are tidally destroyed be- 
fore they can merge are assumed to contribute their stars 
to a "diffuse stellar component" (DSC), which may be as- 
sociated with the stellar halo or the Intra-Cluster Light. 
We have verified that our model reproduces the condi- 
tional multiplicity function of sub-halos over the relevant 
range of host halo masses. Details and tests of our new 
algorithm for the treatment of sub-structure will be pre- 
sented in Maulbetsch et al. (in prep). 



2.2 Cooling 

The rate of gas condensation via atomic cooling is 
computed base d on t he model originally proposed by 
IWhite fc Freii^ (|l99lh . and utilized in various forms in 



virtually all semi-analytic models. Here we use a slightly 
difi'erent implementation of the cooling model than that 
used in SP99 and subsequent papers, which we find is nu- 
merically better behaved. We first compute the "cooling 
time", which is the time required for the gas to radiate 
away all of its energy, assuming that it all starts out at 
the virial temperature: 



^cool 



^fimpkT 



(1) 



Here, /xmp is the mean molecular mass, T is the virial 
temperature TVir ~ 35.9(K'ir/(km/s))^ K, pg{r) is the 
radial density profile of the gas, A(r, Zu) is the tem- 
erature and m ctallicity dependent cooling function 
Sutherland fc Dopita 1993| ) , and Zu is the metallicity 
of the hot halo gas. We assume that the gas density pro- 
file is described by that of a singular isothermal sphere: 
Pai''') ~ '^hot/(47rrvirr'^). Substituting this expression for 
pg{r), we can solve for the cooling radius rcooi, which is 
the radius within which all of the gas can cool within 
a time tcooi. Writing the expression for the mass within 
rcooi, and differentiating, we obtain the rate at which gas 
can cool: 



drrir, 



dt 



1 '■cool 1 
7T"^hot -; 



(2) 



There are various possible choices for the cooling 
time tcooi. Som e early works used th e Hubble time, 
icooi = iff (e.g. iKauffmann et aLlll993l ). In our earlier 
models (e.g. SP99), we used the time since the last halo 
major merger tmrg, defined as a merger in which the halo 
gro ws in mass by a t leas t a fa ctor of two. Here, w e fol- 
low ISpringel et al] (|200lD and ICroton et aP (|2006l ) and 
assume that the cooling time is equal to the halo dy- 
namical time, tcooi ~ tdyn = r^ir/Vvir- Note that be- 
cause in general fdyn < tmrg < tn, and the cooling rate 
dnicooi/dt oc tcooi"^^^, the choice fcooi = tdyn results in 
higher cooling rates than assuming fcooi ~ in , while using 
icooi = imrg produces intermediate results. 

It can occur that rcooi > ^vir, indicating that the 
cooling time is shorter than the dynamical time. In this 
case, we assume that the cooling rate is given by the rate 
at which gas can fall into the halo, which is governed by 
the mass accretion history. 

We note that there are several rather arbitrary 
choices that must be made in any semi-analytic cool- 
ing model — for example, the profile of the hot gas and 
whether it is "reset", and the time to which the cool- 
ing time is compared (see above) — and different groups 
tend to make slightly different choices for these ingredi- 
ents. Changing these ingredients in reasonable ways leads 
to overall variations in the cooling rates (changes in the 
redshift and halo mass dependence tend to be small) of 
at most a factor of two to three. These differences are 
then typically compensated by adjusting the supernova 
feedback and/or AGN feedback parameters. 

We hav e ado pted choices similar to those of 
ICroton et all (|2006l ). in part to facilitate comparison with 
their results, and also because it has been shown that 
this recipe produces good agreement with the cooling 
rates and accumulation of gas in fully 3-D hydrodynamic 
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simulations (|Yoshida et al.ll2002l ) without star formation, 
SN feedback, or chemical enrichment. We have also com- 
pared our results with the cooling rates presented by 
iKeres et al.l (|2005l '). and find good agreement. 

Recently, stud ies based on ID and 3D hydrody- 
nami c simulations (|Birnboim &: Dekelll2003l : iKeres et al.l 
120051 ) have highhghted a distinction between gas which 
is accreted in a "cold flow" mode, in which the gas 
particles are never heated much above ~ 10* K, and 
a "hot flow" mode in which gas is first shock heated 
to close to the virial temperature of the halo, forming 
a quasi-hydrostatic halo, and then cools in a manner 
similar to a classical cooling flow. The possible impor- 
tance of distinguishing between gas flows occuring in the 
regime fcooi < iff vs. icooi > ts (where tg is the free- 
fall time ) has been highlighted many times in the lit- 
erature (Isii5ll977l: lBinnevlll977l : [ Rees fc Ostrikej[l977 : 
IWhite fc Frenklll99ll) . Although other criteria have been 
proposed (see ICroton et al.l l2006l ). we will identify gas 
cooling which occurs in timesteps in which rcooi > f'vir as 
"cold mode" and the reverse (rcooi < r^ii) as "hot mode". 
This distinction will be relevant later, when we begin to 
consider the impact of heating by AGN-driven radio jets. 

As in most semi-analytic models, we assume that all 
new cold gas is accreted by the central galaxy in the halo. 
Because of this, satellite galaxies tend to consume their 
gas and become red, non-starforming, and gas poor. Re- 
alistically, satellite galaxies can probably retain their hot 
gas halos, and thus receive new cold gas, for some time 
after they merge with another halo. This aspect of the 
modelling should be improved; however, for the moment, 
we simply keep this problem in mind, and in some cases 
restrict our analysis to central galaxies. 



the simulations of GOO, reionization occurs fairly late 
(^overlap ~ 8, Zroion = 7). Receut rcsults from the WMAP 
satellite, however , suggest an earlier epoch of reioniza- 
tion, Zrcion ^ 10 (|Spergel et al.l [2OO7I) . We make use of 
the fitting func t ions (B 2) and (B3) from Appendix B of 
iKravtsov etHI (|2004bh to compute the initial fraction of 
baryons that can collapse as a function of halo mass and 
redshift, with Zpv eriap = 11 and gpy eriap — 10. 

As shown bv lSomervilQ (|2002l ) using a similar treat- 
ment of photo-ionization squelching, we find that our 
model reproduces the luminosity function of satellite 
galaxies in the Local Group (Maccio et al. in prep). 



2.4 Disk Sizes 

When gas cools, it is assumed to initially settle into a thin 
exponential disk, supported by its angular momentum. 
We assume that the gas has aquired angular momentum 
before it s collapse, alon g with the dark matter, via tidal 
torques l|Peeblesl [l969l ). Given the halo's concentration 
parameter cnfw, spin parameter A and the fraction of 
baryons in the disk /disk, we can use angular momentum 
conservation arguments to compute the scale radius of 
the exponential disk after collapse. We include the "adi- 
abatic contraction" of the halo due to the gravitational 
force of th e collapsing baryons. O u r approach is base d 
on work bylBlumenthal et al.l l| 19861 ), iFlores et all (|l993l ). 
and'Mo et al." ('1998", MMW98), and is described in detail 
in Somerville ct al. (2008, S08). In S08, we showed that 
this model produces good agreement with the observed 
radial sizes of disks as a function of stellar mass both 
locally and out to 2: '--^ 2. 



2.3 Photo-ionization Squelching 

Photoionization heating may "squelch " or suppress the 
collapse of g as into small mass halos ('Efsta thioulll992l : 
iThoul fc Weinberg 1996. : Quinn et al. 1996). This mecha- 
nism may play an important role in reconciling the (large) 
number of small-mass satellite halos predicted by CDM 
with the o bserved number of satellite galaxies i n the Lo- 
cal Gr oup (|Somervillell2002l : iBenson et al.ll2002l ). iGnedinI 
(|200d . GOO) showed that the fraction of baryons that can 
collapse into halos of a given mass in the presence of a 
photo-ionizing background can be described in terms of 
the "filtering mass" Mf- Halos less massive than Mf 
contain fewer baryons than the universal average. GOO 
parameterized the collapsed baryon fraction as a func- 
tion of redshift and halo mass with the expression: 



/b, coll (2, A'/vir) — 



[1 +0.26MF(z)/Mvirl3 ' 



(3) 



where ft is the universal baryon fraction and Afvir is the 
halo virial mass. 

The filtering mass is a function of redshift, and this 
functi on depends on t h e reion ization history of the uni- 
verse. iKravtsov et"al] (|2004bh provide fitting formulae 
for the filtering mass in the simulations of GOO, pa- 
rameterized according to the redshift at which the first 
-ffii regions begin to overlap (^overlap) and the redshift 
at which most of the medium is reionized (zreion). In 



2.5 Star Formation 

2.5.1 Quiescent Star Formation 

During the 'quiescent' phase of galaxy evolution (i.e. 
in undisturbed disks) we adopt a star formation 
recipe based on the e mpirical Schmidt-Kennicutt law 
(|Kennicuttll 19891 . [ 1998h . The star formation rate density 
(per unit area) is given by: 



Nk 



(4) 



the surface density of cold gas in the disk (in units of 
Mqpc"'^), and Esfr has units of MQyr^^ kpc^^. The 
normalization quoted above is appropriate for a Chabrier 
IMF, and has b een converted from the value given in 
iKennicuttI (|l998l ). which was based on a Salpeter IMF. 

We assume that the gas profile is also an exponential 
disk, with a scale length proportional to the scale-length 



of the stellar disk: 



= Xgas ^disk, where the stellar 



scalelength rdisk is determined as described in tj2.4l We 
adopt Xgas = 1.5, which yields average gas surface densi- 
ties in the ran ge ~ 4 — 60 M^p c"'^, consistent with the 



observations of Kennicufrg (|l998l ) . This value is also con- 
sistent with ob servations of the radia l extent of Hi gas in 
spiral galaxies (|Broeils fc Rhe6lll997l ). 

We further adopt a critical surface density threshold 
Ecrit, and assume that only gas lying at surface densities 
above this value is available for star formation. We can 
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then compute the radius within which the gas density 
exceeds the critical value; 



In 



(5) 



where Eo = mcoid/(27rr-gas). The fraction of the total gas 



supply that is eligible for star formation is then: 



/gas(r < rcrit) = 1 - (1 + rcrit/r-gas)- 



So 



(6) 



and the total star formation rate is: 



EsFR 27rr dr 



1 + 



N Kr cr 



exp(-iVK?'crit/rgas) 



ISchavd (|2004l ') investigated star formation thresholds 
in models of isolated self-gravitating disks embedded in 
dark matter halos, containing metals and dust, and ex- 
posed to a UV background. They found that the gas was 
able to form a cold interstellar phase only above a criti- 
cal surface density threshold of Ecrit ~ 3 — 10 MqPc"'^, 
which is consist ent with observations of SF thresholds in 
spiral galaxies (jMartin fc Kennicutt|[200ll ). We find that 
adopting a value of Ecrit = 6 M0 pc~^ produces good 
agreement with the o bserv ations of global SFR vs. gas 
density of iKennicuttI l|l998h . and also with observed gas 
fractions as a function of stellar mass. 

We allow the normalization of the star formation law 
^Kenn to be adjusted as a free parameter. We find that 
using a value of ^Kenn = 8.33 x 10~^, a factor of two 
lower than the one measured bv IKennicuttI (|l998l ) gives 
good agreement with observed star formation rates and 
gas fractions as a function of stellar mass. We adopt the 
observed value for the slope of the SFR law, Nk ~ 1-4. 

We account for the mass loss from stars (recycled 
gas) using the instantaneous recyling approximation. 
Thus, for an instantaneous star formation rate m, we 
form a mass dm, = {1 — R)mdt of long-lived stars in a 
timestep dt. We adopt a recyc led fraction R = 0.43, ap - 
propriate for a Chabrier IMF l|Bruzual fc CharlotI [20031 '). 



2.5.2 Merger-Driven Starbursts 

As in SPFOl, we parameterize the efficiency of star for- 
mation in a merger-triggered "burst" mode as a function 
of the mass ratio of the merging pair. This is supported 
both by obser vations of star format ion enhancement in 
galaxy pairs (jWoods fc GelleJ 120071') and by nu merical 
simulations of galaxy mergers (ICox et al.l |2003 ). How- 
ever, first, an important question arises: which quantity 
should we use for the mass ratio? Many previous works 
have either used the ratio of the virial masses of the two 
dark matter halos, or else the baryonic masses of the two 
galaxies. Because the ratio of baryons to dark matter can 
vary by several orders of magnitude across halos of dif- 
ferent masses, and moreover is a systematic function of 
halo mass (see i]3.Hl . we find that our results can depend 
quite sensitively on this choice. Moreover, the baryonic 



mass ratio is sensitive to the modelling of supernova and 
AGN feedback. 

When we consider that the simulation results clearly 
indicate that the efficiency of the starburst is mainly de- 
termined by the strength of the torques during the later 
stages of the merger, it is clear that what should be rel- 
evant is the total mass (baryons and dark matter) in the 
central parts of the galaxies. Therefore we define mcoro = 
Muuir < 2rs), i.e., the dark matter mass within twice 
the characteristic NFW scale radius Vs = fvir/cNFW, as- 
suming that the dark matter follows an NFW profile. For 
a Milky Way sized halo (Mvir ~ 2 x IO^^Mq), ~ 27 
kpc and so 2rs corresponds to about 60 kpc, close to the 
scale that we expect to be relevant. We then define the 
mass ratio ^ = (mcoro.i + mbar,i)/(mcore,2 -1- mbar,2), i.e. 
as the ratio of the dark matter "core" plus the total bary- 
onic mass (stars plus cold gas) of the smaller to the larger 
galaxy. 

Now defining eburst as the fraction of the total cold 
gas reservoir in the galaxy that is consumed by the burst, 
we parameterize the burst efficiency via: 



fiburst — fiburst.O /i 



(7) 



This functional form has been shown to describe well 
the scaling of burst efficiency with merger m ass ratio in 
hydr odynamic simulations of galaxy mergers (|Cox et al.l 
I2OO8I . COS). A gain based on COS, we assume that mergers 
with mass ratios below 1:10 do not produce bursts, i.e. 
eburst = for /i < 0.1. 

Numerical studies (|Mihos fc HernquistI Il994l : 
ICox et ahlbOOSh have furthermore shown that the burst 
efficiency in minor mergers (/i ^ 0.25) depends on the 
bulge-to-total ratio of the progenitor galaxies, because 
the presence of a bulge stablizes the galaxy and reduces 
the efficiency of the burst. To refiect the joint dependence 
on merger mass ratio and bulge fraction, we adopt the 
results of COS: 



'O.6I B/r<0.0S5 
7burBt = { 0.74 0.0S5 < B/T < 0.25 
1.02 B/T > 0.25 



(8) 



where B/T is the ratio of the stellar mass in the 
spheroidal component to the total stellar mass (disk plus 
spheroid) in the larger progenitor galaxy at the beginning 
of the merger. 

We have studied the burst efficiency eburst, and 
burst timescale Tburst in equal mass mergers in a large 
suite of numerical simulations containing stellar feedback 
as well as feedback from energy released by accretion onto 
a central black hole (Robertson et al. 2006bi) . We find 
that the burst efficiency can be fit by: 

eburst,o = 0.60(Kir/(km/s))°-''"(l -HgBos)""-''' 

xil + f,r\l + zr-°' (9) 

with a scatter of 4.9 percent, and the burst timescale (as- 
suming a double exponential form for the star formation 
rate peak) is fit by: 

Tburst = 191 Gyr (Kir/(km/s))-'-'*"(l + qEOsf'^ 

il + hr°-'\l + z)-'-'' (10) 
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with a logarithmic scatter of 0.36. Here, V^vir is the virial 
velocity of the progenitor galaxies, qeos is the effective 
equation of state of the gas (see Robertson ct al. 2006b) , 
fg = mcoid / (ficoid + m.star) IS the cold gas fraction in the 
disk, and z is the redshift for which the progenitor disk 
models were constructed. It is important to note that the 
simulations used to obtain these fitting functions span the 
range Kir = 60 - 500 km/s, qsos = 0.25-1, /9=0. 01-0.8, 
and z = 0-6. The results of the fitting formulae should 
be used with caution outside of this range of values for 
the input parameters. 

The parameter qsos can be thought of as parameter- 
izing the multi-phase nature of the Interstellar Medium 
(ISM), such that qEOS ~ corresponds to an isother- 
mal gas, and qsos ~ 1 corresponds to the fully pres- 
surized multiphase ISM. Increasing qsos (or adopting a 
"stiffer" equation of state) increases the dynamical stabil- 
ity of the gas, and suppresses the starburst. Thus larger 
values of qsos give smaller values of eburst.o and larger 
values of rburst (because the burst is more extended). In 
this work we adopt a value corresponding to a stiff equa- 
tion of state, qEOS = 1- 

For the burst efficiency, we can see that the only sig- 
nificant dependence on these parameters is on the equa- 
tion of state qEOS- The burst timescale is more sensitive 
to other parameters (as also found by COS), and in par- 
ticular has quite a strong dependence on Vvir- For our 
adopted fiducial value of qsos ~ 1, the typical value of 
the burst efficiency is eburst ~ 0.8, and the burst timescale 
(exponential decline time) for a Milky Way sized galaxy 
(Kir ~ 130 km/s) is rburst ~ 100 Myr. 

We now parameterize the "burst" mode of star for- 
mation as m« = mburst/Tburst. At the beginning of the 
merger, we allocate a reservoir of "burst fuel" mburat ~ 
Cburst cricoid, where rricoid is the combined cold gas from 
both of the progenitor galaxies. The burst continues until 
this fuel is exhausted, and in the absence of new sources 
of fuel, the burst SFR will decline exponentially, with 
exponential decline time Tburst. However, particularly in 
the early universe, it can frequently happen that a new 
merger occurs while a burst from an earlier merger is 
still going on. In this case, we add the new burst fuel 
to the reservoir, and assign a new burst timescale based 
on the updated galaxy properties. Note that the "quies- 
cent" mode of star formation still goes on as before (the 
burst efficiencies computed from the simulations have the 
quiescent star formation subtracted out). 



2.6 Merger Remnants and Morphology 

2.6.1 Spheroid Formation 

Numerical simulations of mergers of galaxy disks have 
also shown that major mergers fi > 0.25 leave behind 
a spheroidal remnant, while smaller mass ratio minor 
mergers (jj, < 0.25) tend to just thicken the disk, per- 
haps driving minor growth of a spheroid via bar instabil- 
ities. Most previous semi-analytic models have assumed 
a sharp threshold in merger mass ratio (e.g. jj, > /eiup — 
0.25 — 0.3) for determining whether the stars after a 
merger are placed in a "spheroidal" component or not. 
However, in reality there will be a continuum, whereby 



larger mass ratios result in more heating and a transfer 
of more material to a dynamically hot spheroidal compo- 
nent. To represent this continuum, we define the function 



/sph 



/cllip 



(11) 



which determines the fraction of the disk stars that is 
transferred to the "spheroid" or bulge component follow- 
ing a merger (as with bursts, we assume that mergers 
with mass ratio /i < 0.1 have no effect). 

Thus, consider a merger of two galaxies with bulge 
masses Bi and B2, and disk masses Di and D2. The mass 
of the new bulge will be Bnow — B1 + B2 + faph{Di + D2), 
and the mass of the surviving disk will be Dnow = (1 — 
fsph){Di + D2). All new stars formed in the burst mode 
are also deposited in the spheroid component. 

In this paper, we assume that all spheroid growth 
is connected with mergers. That is, we do not consider 
formation of spheroids via disk instabilities. We have ex- 
perimented with including spheroid formation via disk 
instabilities, and find that it has only a minor impact on 
the results presented here. 



2.6.2 Formation of Diffuse Stellar Halos 

There is now considerable observational evidence for spa- 
tially extended stellar components surrou nding bright- 
est g roup and cluster g alaxies (e.g. IZibetti et al.l 
I2OO5I : iGonzalez et all l2005l i. It is now thought that 
these "diffuse stellar halos" (DSH) originated from 
tidally disrupted merging satellit es and /or scattering 
of stars during m a jor mergers (iMurante et al] 2004 • 
iMonaco et"all I2O O6': 'Con rov et al] l2007l : iMurante et al] 
2007': IPurcell et a l. 2007). In our models, the stars from 
all satellites that are deemed to be tidally destroyed be- 
fore they merge (according to the criteria described in 
H2.ip . are deposited in a DSH component. In addition, 
when two galaxies merge, we assume that a fraction 
/scatter of the stars from the satellite may be scattered 
into the DSH (thus the galaxy's mass following a merger 
increases by (1 - /scatter )mstar, sat, where mstar,sat is the 
stellar mass of the merging satellite). 



2.7 Supernova Feedback 

Cold gas may be ejected from the galaxy by winds driven 
by supernova feedback. The rate of reheating of cold gas 
is given by: 



TTlrh 



SJV 



f Kis: 



V2OO km/s 



(12) 



where and Qrh are free parameters (we expect Qrh — 
2 for "energy driven" winds; see e.g. iKauffmann et aP 
I1993I ). We take the circular velocity of the disk Vdisk to 
be equal to the maximum rotation velocity of the DM 
halo, Kiax. 

The heated gas is either trapped within the potential 
well of the dark matter halo, so deposited in the "hot gas" 
reservoir, or is ejected from the halo into the "diffuse" 
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Intergalactic Medium (IGM). The fraction of reheated 
gas that is ejected from the halo is given by: 



Reject ( Vvi 



[1.0 + (Kir/Kjcct) 



acioctl-1 



(13) 



where Oojcct = 6 and 14joct is a free parameter in the 
range ~ 100 -150 km/s. 

We keep track of this ejected gas in a "diffuse 
gas reservoir", which recoUapses into the halo in later 
timesteps and once again become s ava ilable for cool- 
ing. F ollowing [Springel et alJ (j200lh and lOe Lucia et al] 
l|2004l '). we model the rate of reinfall of ejected gas by: 



^^rcinf'all Xreinfall 



'Reject 
^dyn 



(14) 



where Xroinfaii is a free parameter, mojcct is the mass of 
ejected gas in the "diffuse reservoir", and fdyn = ''vir/Vvir 
is the dynamical time of the halo. 

Varying Xreinfaii is degenerate with variations in 
the other supernova feedback parameters, eq^, Qrh, and 
\4ject. For larger values of Xroinfaii, the ejected gas is 
reincorporated in the halo more quickly and tends to 
cool rapidly, so the supernova feedback must be made 
more efficient in order to retain good agreement with the 
abundance of low-mass galaxies. On the other hand, if 
Xroinfaii = (ejected gas is never reaccreted), then the 
baryon fractions in clusters are too low. We have chosen 
to adopt the minimal value of Xroinfaii that allows us to 
fit the cluster baryon fractions and the mass function of 
low-mass galaxies simultaneously. 



2.8 Chemical Evolution 

We track the production of metals using a simple ap- 
proach that is comm only adopt ed in se mi-analytic mod- 
els (see e.g.ISom orvil le fc Prima ck 1999; Cole et al.ll2"o0ol : 
iDe Lucia et"aLH2004l ). In a given timestep, where we cre- 
ate a parcel of new stars dm*, we also create a mass of 
metals AMz ~ ydm,, which we assume to be instanta- 
neously mixed with the cold gas in the disk. The yield 
y is assumed to be constant, and is treated as a free pa- 
rameter in our model. We track the mean metallicity of 
the cold gas Zcoid, and when we create a new parcel of 
stars they are assumed to have the same metallicity as 
the mean metallicity of the cold gas in that timestep. Su- 
pernova feedback ejects metals from the disk, along with 
cold gas. These metals are either mixed with the hot gas 
in the halo, or ejected from the halo into the "diffuse" 
Intergalactic Medium (IGM), in the same proportion as 
the reheated cold gas. The ejected metals in the "diffuse 
gas" reservoir are also reaccreted into the halo in the 
same manner as the gas (see i|2.7|l . 

Throughout this paper, the yield y and all metal- 
licities are given in solar units, which we take to be 
Zq = 0.02. Although this formally represents the total 
metallicity, we note that as we track only the enrichment 
associated with Type II supernovae, our metallicity es- 
timates probably correspond more closely with a-type 
elements. 



2.9 The Growth of Supermassive Black Holes 

We assume that every top-level halo in our merger tree 
contains a seed black hole with mass Mgcod- Typically, 
we assume Msccd — IOOMq, however, we have checked 
that the results presented here are not sensitive to this 
choice for a range of values M^cod ~ 100 - lO^M©. Black 
holes of approximately this mass could be left behind 
as re mnants of massive Pop III stars (e.g. I Abel et"al] 
|2002| ). or could form via direct core-collapse. In our 
models, all "bright mode" accretion onto supermassive 
black holes is triggered by galaxy-galaxy mergers, and 
we assume that this mode of BH accretion is regu- 
lated, and eventually halted, by feedback from the BH 
itself. Our treatment of BH growth and AGN activity is 
closely based on an analysis of a large suite of numeri- 
cal hydro dynamic simulations incl udi ng BH growth and 
feedback (ikobertson et"al] l2006bl l3lal: ICox et al.l l2006bl : 
iHopkins et a l."2007b^ which utilize the m ethodology de- 
veloped in | Di Matteo et al.l (|2005|) and ISpringel et al.l 
(l2005bh . We now briefly summarize the results of those 
simulations and the manner in which we implement them 
in our semi-analytic model. 

In the merger simulations, as the galaxies near their 
final coalescence, the accretion onto the BH rises to 
approximately the Eddington rate. This rapid accre- 
tion continues until the energy being deposited into the 
ISM in the central region of the galaxy is sufficient 
to significantly offset and eventually halt a c cretio n via 
a pressure-drive n outflo w. iDi Matteo et al.] (|2005t ) and 
[Robertson et al.1 l|2006d ) found that the merger simula- 
tions naturally produced black holes and spheroidal rem- 
nants that obeyed the observed BH mass vs. spheroid 
mass relationship. The normalization of the relationship 
depends on the fraction of the AGN's energy that is 
coupled to the ISM, and was chosen to reproduce the 
normalization of the observed r elation. Based on fur - 
ther analysis of these simulations, iHopkins et al.l l|2007bl ) 
suggested that the BH mass is largely determined by 
the depth of the potential well in the cent r al regio ns of 
the galaxy . As shown by iRobertson" et al.l (|2006bh and 
ICox et al.l (2006b"), mergers of progenitor galaxies with 
higher gas fractions suffer more dissipation, and produce 
more compact remnants than those with less gas. There- 
fore, mergers with high gas fractions will produce a rem- 
nant with a deeper potential well and a larger BH mass 
to spheroid mass ratio than gas-poor mergers. This pic- 
ture predicts that there should be a "Black Hole Fun- 
damental Plane" , whereby galaxies with smaller effective 
radius for their mass host larger mass black holes; there 
is observational evidence for the existence of such a BH 
fundamental plane in nearby dormant BH h ost galaxies 
(iMarconi fc Huntll2003l: IHopkins et al.ll2007d ). 

iHopkins et al.1 l|2007bl ) find that the relationship be- 



tween progenitor gas fraction and the final BH mass to 
spheroid stellar mass ratio at the end of the merger ob- 
tained in their simulations can be parameterized as: 

log(MBH/Meph) = -3.27 + 0.36 erf [(/gas - 0.4)/0.28] 

(15) 

with a scatter around this relationship of 0.2-0.3 dex, 
corresponding to the expected range of orbital parame- 
ters. 
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In our semi-analytic model, at the beginning of each 
merger above a critical mass ratio (pcrit ~ 0.1), we 
compute the expected mass of the spheroid that will 
be left behind at the end of the merger, where we as- 
sume that all of the new stars formed in the burst mode 
will end up in the spheroid, along with the 'heated' 
disk stars specified by eqn. [TT] We then use eqn. 1151 
above, to compute mBH, final, the BH mass at the end of 
the merger, based on the initial "effective" gas fraction 

/gas.eflf = ("Icold.l + "lcold,2)/("lbar,l +m.bar,2) (i.e., the 

sum of the cold gas masses in both galaxies, divided by 
the sum of their baryonic masses) . We allow the value of 
JTi-BH, final givcu by Equ. llSl to be scaled by an adjustable 
free parameter /bh, final- 

We assume that the BH in the two progenitor 
galaxies merge rapidly to form a new BH, and that 
mass is conserved in the BH merger. We allow the 
BH to grow at the Eddington rate until it reaches 
a mass A/BH.crit, whereupon it enters the "blowout" 
phase and begins a power-law decline in the accretion 
rate, accordin g to th e family of lightcurves defined by 
iHopkins et aD (|2006bh . From the simulations, MBH,crit ~ 
/BH,crit 1.07 (MBH.finai/lO^Mo) ' , where we introduce 
the adjustable parameter /BH.peak, which determines how 
much of the BH growth occurs in the Eddington-limited 
vs. power-law decline ( "blow-out" ) phases. When the BH 
reaches the mass Mbh. final, "bright mode" accretion is 
switched off. If the pre-existing BH is more massive than 
MBH.crit, it goes straight into the "blowout" mode until it 
reaches Mbh, final. If the pre-existing BH is more massive 
than Mbh, final, the BH does not grow at all, and there is 
no AGN activity. 



2.10 AGN-driven galactic scale winds 

In the numerical merger simulations, the energy be- 
ing released during the rapid growth of the BH also 
drive s powerful galactic-sc ale winds ()Di Matteo et all 
I2OO5I : ISpringel et alll2005ah . We again make use of the 
simulations to parameterize this process in our semi- 
analytic model. We start by equating the momentum as- 
sociated with the radiative energy from the accreting BH 
with the momentum of the outflowing wind: 



1 1 



10 



5. " • • i 



7 8 
log[MBn/Mo] 



Figure 1. Mass ejected in an outflow as a function of the final 
BH mass, from the numerical merger simulations. The dashed 
line shows the scaling predicted by the momentum conserva- 
tion argument (Eqn. I17I I. with twind = 0.5. The dotted line 
shows a simple scaling of ejected mass with BH mass, which 
does not fit the simulation results as well. 
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= M 



outtlow y CSC ) 



(16) 



where twind is the effective coupling efficiency, -Ebh = 
r^ra.dTfia.ccC^ , Moutflow is the mass of the ejected gas, and 
Ksc is the escape velocity of the galaxy. We then obtain 
the following expression for the mass outflow rate due to 
the AGN driven wind: 



Figure 2. Black hole accretion rate as a function of ^ = 

(T7/A23)(m,BH,8), from the observational analysis of Allen et 
al. (2006) (filled circles). The dashed line shows the scaling 
predicted by the NFOO "isothermal cooling flow" model (see 
text), while the dotted line shows a linear flt to the data points. 



dt 



^wind IJrad r- TTZacc ■ 
V^csc 



(17) 



We flnd that this simple formula provides quite a 
good description of the outflow rates in the simulations, 
as shown in Fig. [T] We see that Eqn. [T7] provides a much 
better description of the simulation results than simply 
assuming Moutfl 



2.11 Radio Mode Feedback 

In addition to the rapid growth of BH in the merger- 
fueled, radiatively efficient "bright mode", we assume 
that BH also experience a low-Eddington-ratio, radia- 
tively inefficient mode of growth associated with efficient 
production of radio jets that can heat gas in a quasi- 
hydrostatic hot halo. We base our fiducial model on the 
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assumption that the "radio mode" is fueled by Bondi- 
Hoyle accretion (|Bondilll952l ): 



TTlBondi — 7r(GMBH)^P0Ca ^, 



(18) 



where po = p{ta) is the density of the gas at the accre- 
tion radius ta, Cs is the sound speed of the gas, and we 
have assumed an adiabatic index of 71 — 5/3 for the 
gas. We adopt t h e isot hermal cooling flow solution of 
iNulsen fc FabianI (|2000l . NFOO), in which thermal insta- 
bilities act to maintain the density such that the sound 
crossing time is of order the local cooling time: 

3 iinipkT 



TA 



2p{rA)HT,Zny ^^^^ 

Here, ta = 2GMbh/Cs is the Bondi accretion radius, K is 
a dimensionless constant which depends on the details of 
the flow, kT is the temperature of the gas, and A(r, Zh) 
is the cooling function. Solving for the density po and 
substituting into egn. 1181 we obtain; 



kT 



(- 



Mb 



(20) 



where we have subsumed all constants into the factor 
w^radio. A similar mo del has also been considered by 



IChurazov et~all (|200d ) and lCroton et all (|2006l '). 

We can test the validity of this model using recent 
observations of the central density and temperature of 
hot X -ray emittin g gas in nine nearby elliptical galax- 
ies by lAUen et al.l l|2006l . A06). Deep Chandra observa- 
tions allowed A06 to obtain measurements or reliable 
extrapolations of the gas properties within one order 
of magnitude of the Bondi radius for eight of the sys- 
tems. Each system also has a measured velocity disper- 
sion, which allows an estim a te of t he BH mass using the 
relation of IXremaine et all (|2002l ). In Fig. [3 we com- 
pare the Bondi accretion rates with the quantity C, = 
(T7/A23)(mBH,8), where we define T^ = T/WK, A23 = 
A(r)/(10"^^ergcm^s-^), and mBH,8 = Mbh/IO^^Mq. 
We use the published values of gas density and temper- 
ature and the BH mass estimates from A06, and assume 
that the hot gas has a metallicity of one-third solar. A 
formal fit gives a slope of 1.23 in C,, but we see that the 
NFOO isothermal cooling fiow model is quite consistent 
with the data. In terms of the scaled quantities Ty, A23, 
and mBH,8, Kradio = 2.25 X 10~^ provides the best fit to 
the data. 

In our fiducial semi-analytic model, we assume that 
whenever "hot mode" gas is present in the halo, the cen- 
tral BH accretes at the rate given by eqn. 1181 but we 
allow Kradio to be adjusted as a free parameter. We then 
assume that the energy that effectively couples to and 
heats the hot gas is given by Lhcat = «heat?7radrhradioC^. 
Assuming that all the hot gas is at the virial temperature 
of the halo Tvir, the mass of gas that can be heated per 
unit time is then 



(21) 



^fcr/(Mmp) 



3y2 ' 

4 vir 



using kT/{fimp) = ^V^^ir- The net cooling rate is then the 
usual cooling rate rhcooi minus this heating rate rhheat . If 
the heating rate exceeds the cooling rate, the cooling rate 
is set to zero. 



We apply this heating term only for timesteps in 
which the halos are cooling in the "hot mode" (see i|2.2p . 
That is, if rcooi > '"vir in a given timestep, we assume that 
the gas is not susceptible to the heating by radio jets, so 
it cools at the normal rate. 



3 RESULTS 

3.1 Properties of nearby galaxies 

If we adopt a specific set of values for the cosmological 
parameters, it is relatively straightforward to answer the 
following question: how must dark matter halo mass and 
galaxy mass (or luminosity) be related in order to rec- 
oncile CDM with observations? Numerical N-body simu- 
lations can now accurately predict the multiplicity func- 
tion of dark matter halos and sub-halos (i.e., the num- 
ber density of halos of a given mass), and it is then 
straightforward to adopt a parametric or non-parametric 
model relating halo properties to galaxy properties, and 
to adjust the model to fit the observed stellar mass func- 
tion or luminosity function of galaxies. T his exercise has 
been carried out in terms of luminosity bv lKravtsov et aT 



2004al ) , and in terms of stellar mass bv lWang et al. 



20061 ) and iMoster"et al. (2008, M08). It has been shown 
that if galaxies and halos are related in this way, one 
then also reproduces the observed corr elation functions of 
galaxies as a funct ion of stellar mass l|Wang et al ] |2006l : 
iMoster et all 120081 ). We show the function /star(Mhaio) 
derived by M08 in Fig. O This quantity is defined as 
/star(Mhaio) = JTistar / (/bMhalo) , Or the galaxy's stellar 
mass divided by the universal baryon fraction times the 
halo mass. For central galaxies, Mhaio is the virial mass 
of the halo. For non-central galaxies, Mhaio is the virial 
mass of the halo just before it became subsumed in a 
larger halo. 

Just by comparing the halo mass function with the 
observed stellar mass function (see Fig.[3]right panel), we 
see that in order to reconcile the DM halo mass function 
predicted by CDM with the observed galaxy stellar mass 
function, star formation must not only be inefficient over- 
all (/star ~ 0.2-0.3 at its peak), but the function must be 
a strong function of halo mass. Apparently, the conver- 
sion of baryons into stars is highly inefficient both in small 
mass halos and in large ones, and this efficiency peaks in 
halos with mass ~ 10^^ M0. The interesting question that 
then arises, of course, is which physical processes are re- 
sponsible for shaping this highly variable efficiency, and 
for setting the characteristic halo mass scale - 10^^ Mq? 

The semi-analytic models can give us some insights 
into this question. We note that our adopted halo mass 
resolution (lO^M© for host halos, W^°Mq for sub-halos) 
means that our simulations should be reliable and com- 
plete for galaxies with stellar masses greater than ~ 
IO^Mq. Below this mass, we cannot accurately resolve 
a galaxy's formation history. If we switch off both AGN 
feedback and SN feedbaclu, /star is far too high and 



^ In this model, we fix the metallicity of the hot gas to be 
one-third of solar for purposes of computing the cooling rates. 
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Figure 3. Left panel: Fraction of baryons in the form of stars as a function of halo mass (for central galaxies) or sub-halo 
mass (for satellite galaxies). The solid green lines show the empirical relation (with 1- and 2-cr errors) obtained by Moster et 
al. (2008; see text). Triangles (brown) show the models with no SN or AGN FB; open dots (red) show the model without AGN 
feedback; pentagons (orange) show the Halo Quenching model; solid (blue) squares show the fiducial (isothermal Bondi) model 
with /scatter = 0, and crosses (purple) show /scatter = 0.4. The dashed lines show the sixteen and eighty-fourth percentiles for 
the fiducial model. The diagonal gray line in the bottom left corner shows the stellar mass corresponding to the smallest galaxies 
that we can accurately resolve, ~ IO^Mq. Right panel: Galaxy stellar mass functions for the same models (dotted (brown) no 
SN or AGN FB; short dashed (red) no AGN FB; dot-dashed (orange) Halo Quenching; solid (blue) fiducial isothermal Bondi 
(/scatter = 0); triple dot-dashed (purple) (/scatter = 0.4). Green li nes with error bars show the observed Galaxy Stellar Mass 
functions derived from SDSS by iBell et al.l ll2003bl . light green) and iPanter et al] | |2007| . dark green). The long-dashed gray line 
shows the DM halo mass function with the masses shifted by a factor equal to the universal baryon fraction. 



too flat below 10^^ M© (the mild decline at the low- 
mass end is due to photo-ionization squelching). We ad- 
just the parameters of our model for SN-driven winds 
in order to match the empirical values of /star below 
Afhaio ~ lO^^M©, and find that we require e%fq ~ 1.3, 
a ~ 2, and Viject ~ 120 km/s. Consulting Eqn. 1121 we 
see that this implies that in large galaxies (Vdisk ~ 200 
km/s), the SN-driven mass outflow rate is comparable 
to the star formation rate, and the outflow rate increases 
fairly strongly with decreasing disk circular velocity Vdisk- 
This normalization is in good agreeme nt with the ob- 
servational results of e.g. iMartinI (|l999j). although it is 
unclear that the strong scaling with circular velocity is 
supported by these observations. Also, these winds can 
escape the potential well of the dark matter halo in halos 
with V^ir £j 120 km/s, which is again consistent with the 
observations of iMartig (Il999l ). 

In models with no feedback from AGN, we can see 
from Fig. |3]that /star does turn over at large halo masses: 
this is because large mass objects have formed more re- 
cently, and have had less time to cool. In some of the 
earliest explorations of galaxy formation in the CDM 
paradigm, it was suggested that this cooling time ar- 
gument could explain the chara c teristic mass scale o f 
galaxies l|Blumenthal et al.l [l98i : IWhite fc ReesI Il978l ). 
However, one can see that the turnover occurs at too 
high a mass, and too much gas cools and forms stars 
in large mass halos. This result is obtained not only 
in semi-analytic m o dels by many different groups (e.g. 
iBenson et akl l2003l : ICroton et all l2006l : ICattaneo etahl 



l2006l). but also in numerical hydrodynamic simulations 
(iBalogh et al]|200ll : iBorgani et al]|2006l : ICattaneo eFZI 
l2007l ). 

For comparison. Fig. [3] also shows the predictions 
for /star(Mhaio) in a very simple implementation of the 
concept of AGN heating, what we shall call the "Halo 
Quenching" (HQ) model. In this model, we simply shut 
off cooling flows when the host halo exceeds a mass of 
Mq = 1.3 X IO^^Mq. This model is based on the idea 
that halos around 10^^ M© lie near the transition be- 
tween the "cold flow mode" (gas cooling more rapidly 
than the free-fall time) and the formation of quasi- 
hydrostatic hot halos ("hot flow mode"), and that ra- 
dio jets from SMBH can easily keep gas hot if it is 
in a quasi-hydrostatic hot halo, but not if it is cooling 
in th e cold flow mode (jPekel fc Birnboimlbood : iBinnevI 
|2004| ). This idea has been previously implemented in a 
full semi-analytic model by ICattaneo et al.1 l|2006l ). and 
was found to very successfully reproduce both the lumi- 
nosity function and magnitude dependent color distribu- 
tions of galaxies. We also find that this model reproduces 
/star (A/halo), aud hcuce the galaxy stellar mass function, 
extremely well. 

We can then compare this with the prediction of our 
fiducial model, in which accretion onto a central SMBH 
is modelled assuming Bondi accretion and the isother- 
mal cooling flow model of NFOO (see i|2.1ip . We adjust 
the scaling factor Kradio in order to reproduce /star(Mhaio) 
as well as possible over its whole range. Of course, in this 
model, the scaling of the heating rate as a function of halo 
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mass, and hence the shape of /star(Mhaio), is determined 
by the isothermal Bondi accretion modeL Indeed, we can 
see that our fiducial model (with /scatter = 0) slightly 
overpredicts /star for large mass halos Mhaio ~ 1O^^M0. 
This results in a small excess of large mass galaxies 
(Afstar ^ 10^^ M0) in the predicted stellar mass func- 
tion (Fig. 131 right panel). Note that at halo masses 
Mhaio 10^^, the results for /star are the same as in the 
model without AGN feedback. The radio heating mode 
is ineffective in small mass halos for multiple reasons: 1) 
low mass halos cool mainly in the "cold flow" mode, and 
we have assumed that the "radio mode" is fueled by hot 
gas and that radio jets can only heat gas that is in a 
quasi-hydrostatic hot halo 2) low mass halos tend to host 
disk-dominated galaxies, which do not contain massive 
black holes. 

How concerned should we be about the galaxies in 
high-mass halos being too heavy in our fiducial model? 
The discrepancy amounts to about 0.15-0.2 dex in stel- 
lar masfl. A number of recent observational studies 
have suggested that the luminosities and therefore stel- 
lar masses of the central galaxies in clusters may be 
underestimated by a significant factor (as much as 1.5 
mags) in surveys such as SPSS and 2MASS llLauer et al.l 
I2OO7I : iDesroches eral]|2007l : Ivon der Linden et al.ll2007l ). 
upon which our local galaxy stellar mass function and 
luminosity function estimates are based. Other studies 
have shown that a significant fraction of the stars in 
these galaxies are distributed in a very extended "halo " 
or envelope (iGonzalez et al.1 l2005l : IZibetti et all bOOSl ). 
One possible origin of this extended "diffuse stellar halo" 

i PSH) is stars that are scat tered to large radii in mergers 
Murante et ai]|2004l . |2007| ) . In order to explore this idea, 
we run a model in which a fraction /scatter of the stars in 
merged satellite galaxies is added to such a diffuse com- 
ponent, which is tracked separately from the main stellar 
body of the galaxy. We show the results of such a model 
with /scatter = 0.4 (probably an upper limit on the physi- 
cally plausible value of this parameter) in Fig. [S] and find 
that in this model, the stellar mass function is reproduced 
extremely accurately. 

As discussed in i]2.6.1l in each galaxy, we track sep- 
arately the stars that have survived in an undisturbed 
disk and stars that have been "heated" by mergers to 
form a spheroid. Thus we can assign a crude morphologi- 
cal type based on the ratio of the mass in the "bulge" 
to that in the "disk", B/T. In Fig. g] we show the 
stellar mass functions divided into spheroid-dominated 
and disk-dominated galaxies, compared with observa- 
tional estima tes similar l y divid ed in terms of morpho- 
logical type l|Bell et al.M2003bl ). for the fiducial model 
with /scatter = 0.4. There is a small excess of massive 
disk-dominated galaxies, which may indicate that there 
is still a small degree of overcooling in our most mas- 
sive halos. There is also quite a large excess of low-mass 
spheroid-dominated galaxies, and a deficit of low-mass 



^ Although we do not present any results in terms of luminos- 
ity in this paper, we see a similar discrepancy in the predicted 
luminosity functions in all bands. Therefore we probably can- 
not ascribe the problem to the stellar mass estimates. 



disk-dominated galaxies. These problems persist even if 
we exclude satellite galaxies from our analysis, and can- 
not be eliminated by simply adjusting the parameter /sph 
without ruining the agreement for massive galaxies. 

Another of the important free parameters in our 
model is the normalization of the star formation recipe 
^Kenn- The strougcst constraint on this parameter is the 
ratio of cold gas to stars in galactic disks. Increasing 
^Konn causes gas to be converted into stars more rapidly 
and leads to lower gas fractions. We compare the pre- 
dicted cold gas fractions /gas = mcoid/(mcoid + mstar) as 
a function of stellar mass in our fiducial model with obser- 
vational estimates in Fig. [S] For the models, we consider 
disk-dominated galaxies {B/T < 0.4) which are the cen- 
tral galaxies in their halo (we suspect that the cold gas 
fractions of satellite galaxies may be too low because we 
assume that all cooling gas is accreted onto the central 
galaxy). We comp are with the observational estimates of 
iBell et ai] |20033) for morphologically late-t ype galaxies 
and with galaxies on the blue sequence from iKannappanI 
([2004). The agreement is good at stellar masses greater 
than ~ 10^'^, but the gas fractions are a bit low for lower 
mass galaxies. Note that if we had not adopted a crit- 
ical density in our star formation law (see H2.5p . i.e. if 
we tune Ecrit — > 0, then the gas fractions in low mass 
galaxies come out far too low and we do not reproduce 
the observed trend that gas fractions are higher in low 
mass galaxies. A further check on the cold gas content 
of our galaxies comes from observations of the Hi and 
H2 mass function. We show the prediction of our fiducial 
model compared with these observations in Fig. O and 
find reasonable agreement, but with a hint of a deficit 
of low-gas-mass galaxies, and a small excess on the high- 
gas-mass end. 

A complementary tracer of star formation is the 
heavy elements locked up in stars. We show the stellar- 
mass weighted mean stellar metallicity as a function of 
stellar mass for galaxies in our fiducial model (/scatter = 
0.4) in Fig. |S] Our predictions may be compared with 
observational estima tes based on SDSS spectra from 
iGallazzi etHI ||2005|). I t is w orth noting that the esti- 
mates of lGallazzi et al.l |200^ effectively measure a com- 
bination of a-process elements and Fe, while our mod- 
elling includes only enrichment due to Type II supernova, 
and therefore our "metallicities " correspond more clo sely 
to a-type elements. Also, the IGallazzi et al. I (|2005l ) es- 
timates are effectively luminosity weighted, not stellar 
mass weighted, and may be systematically biased towards 
higher values for supersolar metallicities (see the discus- 
sion m IGallazzi et "aLll2005l l. Considering these potential 
biases, and the relatively crude nature of our chemical 
evolution model, we find fairly good agreement with the 
observed stellar mass vs. metallicity relation. Note that 
although the normalization of this relation can be ad- 
justed by tuning the value of the stellar yield, y, (the re- 
sults shown here adopt y = 1.5 in solar units), the shape 
of this relation is a fairly complex product of various in- 
gredients of the model. For example, the low-mass slope 
is primarily determined by the mass dependent star for- 
mation efficiency caused by our star formation threshold 
and the strongly mass-dependent supernova feedback ef- 
ficiency that we have assumed. The turnover on the high 
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Figure 4. The galaxy stellar mass function divided by morphological type. Solid lines with error bars show observational estimates 
from SDSS (light green) and 2MASS (dark green; Bell et al. 2003) for early type galaxies (left) and late type galaxies (right). Light 
(gray) lines show the observed mass function for galaxies of all morphological types. Dashed lines show the model predictions for 
bulge dominated (B/T > 0.4; left) and disk-dominated {B/T < 0.4; right) galaxies, for the fiducial isothermal Bondi model with 

/scatter = 0.4. 





Figure 5. Left: the cold gas fraction as a function of stellar mass. The (green) squares show observational estimates for mor- 
phologically late-t ype galaxies d erived from the data of iBell et al.l ||20033), and open circles show the observational estimates for 
blue galaxies from iKannappanl | |2004) . The shaded area shows the conditional probability distribution P(/gas|'Tistar) for central 
disk-dominated galaxies predicted by our fiducial (isothermal Bondi) model. The (light blue) solid line shows the median of this 
distribution, and dashed (light blue) lines show the 16 and 84th percentiles. Right: the galactic cold gas mass function. The thick 
(d ark blue) line shows the prediction of our fiducial model. The thick solid ( dark green) curve shows the obs erved Hi mass function 
of lZwaan et al. I Hooi), the da shed (light b l ue) li ne the Hi mass function of lRosenberg &: Schneiderl | |2002|), the dotted (magenta) 
line th e H2 mass function of iKeres et al.l | |2003|) . and the solid (light green) line shows the sum of the iRosenbere fc Schneiderl 
||2002| ) Hi mass function and the H2 mass function. 
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Figure 6. Left: Stellar mass vs. stellar metallicity. Grey shading shows the conditional probability P(Z»|m») for our fiducial 
model (/scatter = 0.4), and the (light blue) solid and dashed lines show th e 50th, 16th and 84th percentiles. The (green) dots 
and dashed lines show the observational estimates from ICallazzi et al. Right: Stellar mass vs. stellar mass-weighted mean 

stellar a ge. Shading, lines, an d symbols are as in the left panel. Open square symbols show the predictions of the semi-analytic 
model of lCroton et al.l ||2006| ') for comparison. 



mass end is shaped by the quenching of star formation 
by AGN feedback and gas-poor mergers. 

In Fig. |6] we also show our model predictions for 
the stellar mass weighted mean stellar age of galaxies 
as a function of stell ar mass, compared w ith the obser- 
vational estimates of iGallazzi et al.1 (|2005l ). The models 
predict a weak trend of older ages in more massive galax- 
ies, with the ages of massive galaxies slightly older than 
the observational estimates. However, the predicted trend 
in our m odels is much weake r than the observed trend 
found bv lGallazzi et al.l l|2005l ). and low-mass galaxies in 
our models are mu ch older than the observations indi- 
cate. ICroton et al. showed that models without 
AGN feedback predicted that massive galaxies have ages 
as young as low-mass galaxies, and that introducing ra- 
dio mode AGN feedback produced an age-mass trend 
with the correct sense (more massive galaxies are older). 
However, they did not compare directly with observa- 
tional estimates. We reproduce the predict ions from the 
AGN feedback model shown in Figure 10 of lCroton et al.l 
(|2006l ). and see that their results are very similar to ours 
(in fact low-mass galaxies are slightly older in their mod- 
els than in ours) . We should keep in mind that if a galaxy 
has a significant-by-mass older stellar population with a 
small "frosting" of young stars, the ages derived from 
stellar absorpt i on lin es (mainly Balmer lines) as in the 
iGallazzi et al] l|2005l ) approach will be biased towards 
young ages. This discrepancy is worth examining in more 
detail, but this is beyond the scope of this paper. 

The relationship between galaxy mass and BH mass 
is clearly a key result that our model should reproduce. 
Recall that in our model, this relationship is set by the 
depth of the potential well of the galaxy at the time when 
the BH forms, which in turn is determined by the gas 




Figure 7. Predicted relationship between bulge mass and 
black hole mass (grey shading indicates the conditional prob- 
ability P(mbhl™bulge); light blue solid and dashed line shows 
the median and 16th and 84th percenti l es) co mpared with 
the observed relation from iHaring fc Rb] ||2004| . green lines). 
Symbols show the m easurements for individual galaxies from 
iHaring fc Rb] ll2004h . 



fraction of the progenitor galaxies of the last merger (see 
ij2.9p . More gas-rich progenitors suffer more dissipation 
when they merge, and produce more compact remnants 
with deeper potential wells. A deeper potential well re- 
quires more energy, and therefore a more massive BH 
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in order to halt further accretion and growth. Although 
we have seen that the predicted gas fractions of disks at 
the present day agree reasonably well with observations, 
the gas fractions of the progenitors of black hole hosts 
depend on many factors, such as the masses of those pro- 
genitors at the time when the BH is formed, the epoch of 
formation of BH of a given mass, and the details of the 
star formation and feedback modelling. It is therefore a 
non-trivial success of our model that we reproduce the 
observed slope and scatter of the MsH-Maph (black hole 
mass vs. spheroid mass) relationship, as seen in Fig. [7] It 
is interesting that our model predict s a smal l upward cur- 
vature at the high-mass end, which IWvithd I2OO6. ') argue 
is present in the observed relation. Our predicted rela- 
tion also has a somewhat flatte r slope at low BH ma sses 
than the extrapolation of the iHaring fc Rixl l|2004 ) re- 
sults; however, there are currently very few robust BH 
mass estimates at such low masses. 




3.2 Group and Cluster Properties 

We have so far focussed on the properties of individual 
galaxies. We now consider predictions for a few properties 
of groups and clusters. We have selected these quantities 
because they help to constrain some of the free param- 
eters or uncertain ingredients in our models. In Fig. [S] 
we show the hot gas fraction (/hot = ?rihot/Mvir), i.e., the 
mass of hot gas contained in the dark matter halo divided 
by the total virial mass of the halo. The hot gas fraction 
in our models has a sharp "step" at Mvir ~ lO^'^M©, be- 
cause of the rapid transition between halos in which the 
gas cools rapidly compared with the dynamical time, so 
there is typically little hot gas present in the halo, and ha- 
los in which the cooling time is longer compared with the 
dynamical time, so the halo can build up a reservoir of hot 
gas. In our fiducial model, this hot gas is then maintained 
by AGN "radio mode" heating. Our results are in reason- 
able agreement with the hot gas fractions in clusters es- 
timated from ob servations of their X-ray emitting gas by 
IVikhlinin et al.l ((2OO6 ) . These observations are somewhat 
uncertain, because the X-ray emission typically cannot be 
detected all the way out to the cluster virial radius, so it 
must be extrapolated. However, these observations pro- 
vide an important constraint on the modelling of re-infall 
of gas that has been ejected by supernovae (see i|2.7|l . If 
we do not allow this ejected gas to be reaccreted at all, 
then the baryon fractions in clusters are predicted to be 
significantly smaller than the universal value, in confiict 
with observations. 

In Fig. [S] we also show the predicted metallicity of 
the hot gas in halos. In our model, the hot gas is en- 
riched by the ejection of metals from the cold gas in 
galactic disks by supernova-driven winds. We tuned the 
chemical yield y to reproduce the metallicities of stars 
in galaxies, so the metallicity of the hot cluster gas is 
a cross-check on our chemical evolution and supernova 
feedback modelling. We find that the hot gas in cluster- 
mass halos is enriched to about 0.25 of the solar value, 
and is nearly constant above about Afvir ~ 10^"^ Mq. This 
is close t o the value of ~ Q. ?>Z(:) measured for hot gas in 
clusters l|Arnaud et al.|[l993 ) . These measurements of hot 



Figure 9. Mass contained in a "diffuse stellar halo" (DSH) 
relative to the mass of the central galaxy plus the DSH mass. 
The shading and (light blue) curves show the predictions of 
our fiducial model; the green bo x shows the approxima te locus 
of observational estimates from lGonzalez et al.l ll2005h . 



gas metallicity are primarily sensitive to Iron, while as we 
discussed earlier in this section, our chemical evolution 
modelling traces only the metals produced by Type II 
supernovae, so we do not expect perfect agreement. Also, 
some additional metals may be dr iven out of the ga laxies 
by strong shocks during mergers (jCox et al ] |2006al ). 

We have argued that a significant fraction of stars 
may be scattered into a "diffuse stellar halo" (DSH) by 
mergers. It is important to check whether the predicted 
mass of stars in these DSH is in agreement with direct 
observational measurements. In Fig. [9] we show the mass 
of the DSH divided by the total mass of the DSH plus the 
main galaxy (/dsh = mDSH/(mDSH + mBCc)), as a func- 
tion of the virial mass of the halo. The model predictions 
are c onsistent with th e rang e of observational estimates 
from iGonzalez et al.l (|2005f ). which we hav e ado pted 
from the results presented in ICon rov et al .l ll200'if). In 
agre ement with pr e vious studies by Monaco et al. ( 20061 ) 
and IConrov et al] l|2007h . we find that the model with 
/scatter = 0.4 is able to reproduce the observed stellar 
mass function of galaxies at z ~ as well as the fraction 
of stellar mass in the DSH component 0. 



^ We note that th e published value of /dsh — 0.33 from 
IZibetti et al. is significantly lower than the values we 

have adopted here. While several factors may play a role, such 
as the different extent of the photometry and s ample selection, 
much of the discrepancy with the results of iGonzalez et al.l 
1 I2OO5I) is apparently due to the details of the way the BCG and 
DSH (or intracluster light) are defined (A . Gonzalez and S. Zi- 
betti, priv. comm.; see also lZibett 3 1I2OO8I '). Section 5.1). When 
Zibetti adopts the same decomposition method as G onzalez, 
he fin ds much more consistent values /dsh — 0.67 llZibettil 
I2OO8I) . 
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Figure 8. Left: the fraction of lialo baryons in the form of hot gas as a function of halo virial mass. The gray shaded area shows 
the conditional probability distribution P(/hot I A/vir) for our fiducial model, with the 16, 50, and 84th percentiles shown with (light 
blue) curves. The sharp discontinuity at M^ir — 10^^ Mq represents the transition from rapid co oling (cold flows ) to th e formation 
of a hot halo. The (green) solid circles show the observational estimates of hot gas fraction from lVikhlinin et al.l {ioO^. Right: the 
metallicity of hot cluster gas as a function of halo m ass in our fiducial m odel. The green dashed line shows the observed value in 
clusters, approximately one-third of the Solar value llArnaud et al.lll992l') . 



3.3 Radio Mode Heating 

Several other groups have implemented heating by radio 
jets from AGN in semi-analytic models, and shown that 
in this way they can so lve the overcooling problem and 



:o olmg pr 

other related proble ms ( Croton et al. I I2OO6I : iBower et al.l 
I2OO6I : [Monaco ct al. 2007 ). However, these works have not 



addressed whether the amount of energy required or the 
scalings as a function of halo mass adopted in these mod- 
els are consistent with constraints from observations of 
radio galaxies and cooling flow clusters. We turn now to 
this question. Fig. 1 101 shows the predicted cooling rate as 
a function of halo mass in a model without AGN heat- 
ing. The shaded area shows the cooling rates predicted 
by the full semi-analytic model, while the smooth lines 
show the cooling rate given by Eqn. [S] assuming that 
^■hot ~ /tAfvir and Zhot/-^© = 0.33 (we refer to this 
as the "static halo" cooling model). The lower of these 
lines is for redshift z — 0, and the the higher is for z — 2. 
The divide between "cold mode" and "hot mode" halos at 
~ IO^^Mq is indicated by the color of the lines (with blue 
indicating cold mode, red indicating hot mode). Clearly, 
the static halo model prediction can differ from the cool- 
ing rate in the full SAM because, as we saw in Fig. [8] halos 
of a given mass have a range of values of hot gas fraction 
and metallicity due to their different formation histories. 
Here we see that, in the absence of AGN heating, clus- 
ter mass halos would be expected to have cooling flows of 
hundreds up to one thousand solar masses per year, which 
we know to be in conflict with X-ray observations. In the 
middle panel we show the rate at which gas is heated by 
the "radio jets" in our fiducial (isothermal Bondi) model. 
Note that although we show non-zero heating rates below 



~ IO^^Mq, actually, most of these halos are cooling in the 
"cold flow" mode and so their cold gas accretion rates are 
unaffected by the AGN heating. We note that the heat- 
ing and cooling rates cross near the "magic" halo mass of 
~ IO^'^Mq, and that the heating rate is a steeper function 
of halo mass than the cooling rate at large masses, so that 
there is a lot of "excess" energy being deposited in the hot 
gas. At the moment, in our simple modelling, this excess 
energy is not accounted for. The final panel in this plot 
shows the net cooling rate including the AGN heating. 
Cooling fiows are quenched entirely in halos more massive 
than Afvir ~ lO^'^Mo. However, this is not a sharp cutoff. 
There is a transition region between 10^^ < Mvir lO" 
where some halos have had their cooling flows quenched 
and some have only had them reduced. 

We can compare the heating rates needed in our 
model in order to solve the overcooling problem and the 
galaxy mass problem with observations of the hot bub- 
bles associated with radio jets, seen in the X-ray gas in 
groups and clusters. By estimating the amount of energy 
required to inflate the bubbles, these systems can be used 
as "calorimeters" , gi ving an estimate o f the power being 
injected by the jets. lAllen et all (I2OO6I ) find a tight cor- 
relation between the Bondi accretion rate, and the jet 
power. For systems that also have a black hole mass esti- 
mate or its proxy, for example from a measured velocity 
dispersion or bulge mass, we can then assess the fraction 
of the black hole's rest mass that is being extracted as 
kinetic energy that can heat the gas. In Fig. 1111 we show 
observational estimates of the rate of energy injection, 
or jet power, as a function of black hole mass, from ob- 
servations of elliptical galaxi es wi th associated hot ga s 
bubbles by I Allen et al] l|2006l ) and lRaffertv et al.l (|2006l ). 
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Figure 11. Left: Rate of energy inp ut by the "radio m ode" (jet power) as a function of BH mass. Large, solid green circles 
show the observational es timates from lXllen et all fcOOSh . Open and skeletal triangles (orange) show the observational estimates 
from lRaffertv et alj ||2006| ), for systems at redshift z < 0.05 and z > 0.05, respectively. The thick (dark green) solid line shows the 
time-averaged heating rate derived from observations bv lBest et al.l ll200(j) . The gray shaded area shows the conditional probability 
distribution P(Pjet I^^bh)! and the light blue curves show the median and 16 and 84th percentiles of this distribution. The (red) 
dashed line shows the heating rate that would result if all BH accreted at a fixed fraction {fn = 3.5 X 10^^ is shown here) of 
their Eddington rate, and the (purple) dot-dashed line shows Pj^t cx; ffiglj. Right panel: the jet power divided by the black hole 

1/2 

mass (specific jet power) as a function of halo virial velocity (oc T^/^ ). The gray shaded region shows the conditional probability 
distribution, and light blue lines show the 2.2, 16, 50, 84, and 98th percentiles (from bottom to top), for our fiducial isothermal 
Bondi model. The dot-dashed (orange) line shows the scaling for the fiducial model of ICroton et al.l ||2006| ). 



The iRafFcrtv et al.' ('2006') estimates of jet power overlap 
those of Allen ct al. (2006), but extend to considerably 
higher values for a given BH mass. This may be because 
the RafFerty et al. sample, which extends to higher red- 
shift, con tains more massive clusters than the very nearby 
sample of lAUen et all (|2006l ). These massive clusters may 
have higher gas densities at the relevant radii, thus al- 
lowing more efficient coupling of the radio jet with the 
ICM (S. Allen, priv. comm.). We also show the time av- 
er aged heating rate as a function of BH mass estimated 
by Best et al. (2006), from observations of the radio-loud 
fraction of SDSS galaxies. The observed scaling of jet 
power with BH mass is a bit steeper than it would be if 
the jet power were proportional to the Eddington lumi- 
nosity: the jet power scales approximately as Mbh^'^, or 
as LEddMBH° ^ 

These results may be compared directly with the "jet 
power" as a function of BH mass incorporated in our fidu- 
cial model. Recall that in our model, we assumed that the 
central density and temperature of the gas was set by the 
isothermal cooling flow model of NFOO (see tj2.1ip and the 
accretion rate onto the BH was then set by the Bondi ac- 
cretion rate. We allowed an overall scaling factor Kradio, 
which we adjusted to the minimum value that produced 
a good fit to the empirical constraint on galaxy stellar 
mass as a function of host halo mass discussed in i|3.1l We 
find that the resulting accretion rates are about an order 
of magnitude high er than the Bondi rates estimated by 
I Allen et aU (|2006l ). and the jet power at a given BH mass 



is also higher than the lAllen et al.l (|2006l ) estimates. How- 
ever, the jet powers are consis tent with the higher val- 
ues in the lRaffertv et ahl l|2006l ) sample. It is also impor- 
tant to remember that these observational estimates are 
lower limits. They include only the energy associated with 
inflating the bubbles, while significant energy can also 
be dissi pated through sound wa v es, viscosity, and we ak 
shocks (iMcNamara et al.l l2005l: iNulsen et all l2005allbl: 



iFabian et all l2006l : iForman et al.l |2007| ). iBinnev et al.l 
( 20071 ) analyzed 3D adaptive grid simulations of heat- 
ing of cooling flows, and found that the bubbles reflected 
only ~ 10 percent of the total injected energy. 



However, it is also possible that we are overestimat- 
ing the energy required, because we are insisting that 
AGN heating does the whole job, while as we have dis- 
cussed, there may be several other processes that help 
to reduce the efficiency of cooling in group and cluster- 
mass halos. Furthermore, the semi-analytic cooling model 
is probably only accurate to a factor of two to three at 
best, and may be overestimating the cooling rates. We 
intend to test and better calibrate our models by compar- 
ing with the results of numerical simulations, but this is 
not entirely straightforward, because supernova feedback 
probably plays an important role in altering the equation 
of state and metallicity of the gas, which affects the cool- 
ing rates. Given these uncertainties, we conclude that the 
heating due to AGN predicted by our simple model is not 
only very successful at solving the overcooling problem, it 
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Figure 10. Top: cooling rate as function of halo mass, for 
model with no AGN feedback. Grey shading indicates the con- 
ditional probability Pirn^aoA^h)- Light blue lines indicate the 
2.2, 16, 50, 84, and 98th percentiles. The smooth red and blue 
solid and dashed lines indicate the expected cooling rates in 
a "static halo" model (see text) for z = and 2 = 2, re- 
spectively; the dark blue section of these lines indicates the 
approximate regime for "cold mode" infall, and the red lines 
for "hot mode" (see text). Middle: heating rate by "radio jets" 
in our fiducial (isothermal Bondi) model. Bottom: net cooling 
rate in fiducial model with radio mode feedback. Halos with 
cooling rates below 10~*'^MQyr~^ are plotted in the bottom- 
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is also reasonably consistent with the direct observational 
constraints. 

Fig. [TT] (right panel) also shows the predicted "spe- 
cific jet energy" (jet energy divided by BH mass) as a 
function of the halo virial velocity. One can see that the 
isothermal Bondi model predicts a fairly strong depen- 
dence of jet power on halo virial velocity (or temper- 
ature), and in fact this scaling is important in fitting 
the shape of the /star(Afhaio) function at intermediate 
masses — we find that if we assume that the accretion 
rate (and hence the jet power) is just a function of BH 
mass, we obtain qualitatively similar results, but we do 
not get as good a match to the /star(Afhaio) function 
and hence the galaxy stellar mass function. For com- 
parison, we also show the empi rical scaling adopte d in 
the fiducial Munich SAM ( e.g. [Croton et al.ll2006l ). As 
also noted bv lCroton et al. I i|2oogr their adopted scaling 
Wacc/wBH oc Vjir is Very similar to that predicted by the 
isothermal Bondi model, and it yields very similar results 
when we adopt it in our SAM. 



3.4 Star formation Quenching 

As we noted in the introduction, there are two puzzles 
in galaxy formation that we wish to address in this pa- 
per. One is that real galaxies do not grow as massive as 
we would have predicted in the absence of AGN feed- 
back. The other is that star formation in most massive 
galaxies has been quenched, while most low-mass galax- 
ies continue to form stars. These two problems seem very 
likely to be interconnected, but it is not obvious that a 
specific model which solves one problem will necessarily 
solve the other. In order to assess the quenching prob- 
lem, it is common practice to compare model predictions 
with observed optical or optical-NIR color-magnitude 
distributions. However, these kinds of predictions are 
quite sensitive to metallicity and dust, which compli- 
cates the interpretation. Instead, we make use of the 
physical properties, specific star formation rate (SSFR 
= ?Tistar/m-star) and Stellar mass, derived from GALEX 
UV photometry plus SPSS five-band optical photome- 
trv (ISalim et aLll2007l:I Schimin ovich et al.ll2007l ) . lYi et al.l 
(120051 ) and iKavirai et al. (2003) have shown that the 
NUV-optical colors are a highly effective way to probe 
small amounts of recent star formation in galaxies with 
an underlying old stellar population. 

Fig. ll2l shows the conditional probability distribution 
of SSFR as a function of stellar mass P(SSFR|mstar). The 
top left panel shows the observed distribution derived 
from the GALEX-I-SDSS data by ISchiminovich et al.l 
(|2007l ). The star forming sequence, quenched sequence, 
and the dividing line (sometimes calle d the "gre e n val - 
ley") derived from GALEX+SDSS by Salim et al.l (|2007l ) 
are also shown, and are repeated on every panel. It is 
important to remember that the GALEX-SDSS survey 
is incomplete in the bottom left region of the plot. We 
also show the same distribution, as predicted by the no 
AGN FB model (top right), the Halo Quenching model 
(bottom left), and the fiducial isothermal Bondi model 
(bottom right). The star formation rates shown for the 
models have been averaged over the past 10* years. 




Figure 12. Specific star formation rate (star formation rate divided by stellar mass) vs. stellar mass. Gray shading and contours 
indicate the conditional probability P(SSFRjmstar)- The diagonal (dark blue) solid line in the upper left part of the plot and the 
(red) line in the lower right part of the plot indicate the "star forming sequence" and "quenched sequence" from the observational 
results of ISalim et al.l |2007.). The middle, dashed (green) line indicates the dividing line between the star forming or active 
galaxies and quenched galaxies (sometimes called the "green valley"). These active, valley, and quenched sequences based on the 
observed distributions from GALEX are repeated on all four panels. Top left: observed SSFR vs. mass distribution from GALEX 
ijSchiminovich et al.ll2007h . Top right: predicted distribution from the model with no AGN feedback. Bottom left: predicted 
distribution from the HQ model. Bottom right: predicted distribution from the fiducial (isothermal Bondi) model. 



All of the models shown have the same problem with 
low-mass galaxies. The star forming sequence is nearly 
fiat, rather than being tilted such that less massive galax- 
ies have higher SSFR, as in the observations, and the spe- 
cific star formation rates are too low. Though we have 
tried extensive experiments with parameter variation, we 
have not succeeded in solving this problem. The problem 
also seems to be quite robust to the star formation recipe 
that we adopt. If we remove the SF threshold, thereby 
effectively increasing the star formation efficiency in low- 
mass galaxies, then the galaxies consume more gas and 
have lower gas fractions at the present day. Their SSFR 
are still low, because they have very little fuel for star 
formation. If we increase the star formation threshold. 



which makes star formation even more inefficient in low- 
mass galaxies, the galaxies have higher gas fractions at 
the present day, but most of that gas is not allowed to 
make stars, so the SSFR are again almost the same as 
before. 

In the model without AGN feedback, nearly all mas- 
sive galaxies have high SSFR and would also have blue 
colors, in drastic conflict with the observations. In the 
HQ model, which did as well or better than our flducial 
model at matching /star(A^haio) and the galaxy stellar 
mass function, essentially all massive galaxies are com- 
pletely quenched. This is not surprising, as we know that 
there is a fairly tight relationship between stellar mass 
and halo mass, so a sharp cutoff in halo mass produces 



© 0000 RAS, MNRAS 000, 000-000 



Co-evolution of Galaxies, Black Holes, and AGN 23 



10" 



-5 10-' 

E 

aT 10"* 
u. 
I/) 



io-» 
io-« 



9.18 


V/ \ 






\\ 
\ \ 



■13 -12 -11 -10 -9 
log SFR/M^ (yr-'l 




■13 -12 -11 -10 -9 
log SFR/M,, [yr-] 



■ 13 -12 -11 -10 -9 
log SFR/M,, (yr-J 




■13 -12 -11 -10 -9 
log SFR/M^ (yr-'l 



■13 -12 -11 -10 -9 
log SFR/M., [yr-] 



10" 



1 10-' 
E 

(£ lO"* 

u. 

(/) 

tn 



10-* 

io-« 



10.72 




i 1 
i 


\V 
\\ 

A 

A 

1 V- 



■13 -12 -11 -10 -9 
log SFR/M„, (yr -J 





10 • 


1 1.03 






10-* 






E 








of 


10-* 




>^^^\ 
- - "V* 


u. 
ID 








io-» 


/ / ,y 


Vi 




io-« 


'' 





■13 -12 -11 -10 -9 
log SFR/M^ (yr-'l 




■ 13 -12 -11 -10 -9 
log SFR/M,, (yr-'l 



10- 



-\ 10-' 
E 

0? 10"* 
u. 
in 
in 

• 10-* 



10- 



1 1.65 










!\ 




'. i 








1 
1 



■ 13 -12 -11 -10 -9 
log SFR/M„, (yr-'J 



Figure 13. Distribution of specific star formatio n rates (SSFR) in stellar mass bins (as indicated on the panels). Solid (green) 
lines show the observational results from GALEX ijSchiminovich et aI.ll20o"^ . dot-dashed (red) lines show the no AGN FB model, 
and dashed (purple) lines show the fiducial (isothermal Bondi) model. Vertical dotted lines show the rough location of the division 
between the "active" and "quenched" populations. 



a fairly sharp cutoff in stellar mass. This model cannot 
account for the population of massive galaxies with small 
but detectable amounts of recent st ar for mation, seen in 
the GALEX observations (Yi et al. |2005| argue that the 
NUV hght in these galaxies is indeed due to star forma- 
tion and not a UV upturn or AGN). However, our fiducial 
isothermal Bondi model does produce such a population 
of massive galaxies whose star formation has been sub- 
stantially, but not completely, quenched. The SSFR vs. 
TTistar distribution predicted by this model looks qual- 
itatively quite similar to the observations (for massive 
galaxies) . 

We analyze the distribution of SSFR vs. stellar mass 
in more detail in Fig. 1131 Here, we show histograms of 
SSFR in stellar mass bins, for the observations and two 
of the models: the no AGN FB model, and the fiducial 



model. From this comparison we can see that our fiducial 
model is not producing quite enough massive, actively 
star-forming galaxies (mstar ^ 10.7). It seems that star 
formation is actually being quenched a bit too effectively 
in massive galaxies. However, overall the agreement is 
quite good. 

3.5 Global formation histories of stars, cold 
gas, and metals 

In this paper, we have focussed mainly on predictions 
of the present-day (z ~ 0) properties of galaxies. We 
plan to explore the predictions for the properties of high- 
redshift galaxies in detail in future papers. However, in 
this section we present predictions for the global histo- 
ries of several important (and observationally accessible) 
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Figure 14. In all panels, dotted (red), dot-dashed (orange), and dashed (blue) lines show predictions from the C-ACDM no AGN 
FB model, the HQ model, and the fiducial (isothermal Bondi) model, respectively. The triple-dot dashed (purple) line shows the 
WMAP3 model. Top left: Star formation rate density as a function of redshift. The upper set of thicker lines shows the total SFR 
in the mod e ls, an d the lower set of thin lines shows the SFR due to bursts. Symbols show the compilation of observational results 
of [Hopkins (2004), converted to a Chabrier IMF. The thick solid (gray) curve is the fit to the observational com pilation presented 
bylHopkins fc Bcacom: (.2006, ) . The thin (gray) solid broken curve shows observational estimates from G ALEX tiSchiminovich et al.l 
2OO5I). Top right: The int egrated globa l stell ar mass density as a functi on of redshift. Symbo ls show observational estimates from 
Bell et al.l ll2003bl . 2 ~ O^. lBorch et al.l 1 I2OO6I. z ~ 0.2-1; diamonds), and lFontana et al.l pOOd _ z ~ 0.2-3.5; filled circles). The thick 
(gray) curve shows the fit to the observational compilation presented in IWilki'ns'er^.l l|2008l ). Bottom left: The mass density of 
cold gas in units of the critic al density Ocoid = Pcold/Pcrit- Symbols show the observational e stimates at 2 ~ from the blind Hj 
survey of lZwaan et al.l ll2005h . and from Damped Lyman-a systems (crosses) at high redshift l lProchaska et ani2005h . Small open 
squares also include the contribution from lower column density absorption systems (Lyman-limit systems). Bottom right: The 
mass- weighted average metallicity of stars (.^star ) = Pz /Pstar a s a fun ction of redshift. The diamond symbol at 2 ~ 0.1 shows the 
observational estimate from SDSS galaxies from lGallazzi et al. I ll2007t) . 



components of galaxies: the star formation rate, stars, 
cold gas, metals, and black holes. 



3.5.1 Dependence on Cosmology: the WMAP3 model 

In this section we introduce a new model, which has the 
same recipes for galaxy and BH formation as our fidu- 
cial isothermal Bondi model, but adopts the cosmologi- 



cal parameters (see Table [T| f rom the 3-year analy sis of 
the WMAP data reported in ISpergel et all (|2007l ). We 
will refer to this as the WMAP3 model Q. For our pur- 

* While this paper was in the refeering process, new results for 
the cosmological parameters derived from the five year WMAP 
data, combined with distance estimates from Type la SNae 
and Baryon Acoustic Oscillations, were posted on astro-ph 
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poses, the most important differences between the Con- 
cordance ACDM (C-ACDM) model that we have been 
using so far and the WMAP3 model is that WMAP3 
has a lower value of ag, the normalization of the primor- 
dial power spectrum, and the primordial power spectrum 
also has a "tilt" (n^ = 0.96) while C-ACDM has a scale- 
free initial power spectrum {us = 1). This results in less 
power on small scales, and hence later structure forma- 
tion in WMAP3 relative to C-ACDM. In order to repro- 
duce the 2 = observations as before, we re-tuned the 
star formation efficiency, radio mode heating efficiency, 
and scattering parameters (we used AKenn = 1-67 x 10"*, 
= 6.0 X 10"^, and /scatter = 0.2), but left the other 
parameters the same. After this re-tuning, the WMAP3 
model produces nearly indistinguishable results from the 
C-ACDM model for all of the quantit ies that we have 
shown so far (see also lWang et al.ll2007l ). 



3.5.2 The Global Star Formation and Mass Assembly 
History 

Fig. [TJ] (top left) shows the global star formation rate 
density of all galaxies predicted in the three C-ACDM 
models: the model with no AGN FB, the Halo Quench- 
ing model, and the fiducial isothermal Bondi model. A 
compilation of observational estimates is also shown. The 
no AGN FB model overpredicts the amount of star for- 
mation at low redshift, and the SFR does not decline as 
rapidly as the data indicate, while both the HQ model 
and the isothermal Bondi models produce very good 
agreement with the observations at z ^ 2. We can see 
that the effect of the AGN heating starts to become sig- 
nificant only ai z ^ 4. At higher redshifts, the AGN heat- 
ing has little or no impact on the global SFR, because 
most star formation is taking place in relatively small 
mass halos, which are not affected by the "radio mode" 
heating, as we have already discussed. It is also inter- 
esting that the simple HQ model produces such similar 
results to the fiducial isothermal Bondi model, even over 
a large range in redshift. At 2 ^ 4, the C-ACDM models 
predict m ore star formation th an the observational com- 
pilation of lHopkins fc Beacomi (2006) , but agree with the 
higher estimates in the literature (e.g. ISteidel et al.ll 19991 : 
iGiavalisco et al.l 120041 '). The WMAP3 model predicts a 
much more rapidly declining SFR at 2 > 2, somewhat 
lower than the observations but well within the obser- 
vational errors. 

The lower set of lines shows the star formation con- 
tributed by merger-triggered bursts. The contribution 
due to bursts is much lower than in our previous mod- 
els (SPFOl), for several reasons. (1) Our new treatment 
of dynamical friction and tidal stripping and destruction 
means that many low-mass satellites take longer than a 
Hubble time to merge, or are tidally destroyed before 
they can merge. Therefore the number of minor mergers 



is lower. Also, we do not include satellite-satellite merg- 
ers here, which were included in our previous models. (2) 
Our newly calibrated burst efficiencies in minor mergers, 
which are based on a greatly expanded and improved set 
of hydrodynamic simulations, are lower than in our previ- 
ous models (3) Our Kennicutt star formation law causes 
the star formation efficiencies in quiescent disks to in- 
crease with increasing redshift. As we already showed in 
SPFOl, this leaves less cold gas fuel for bursts, leading to 
a decreased contribution to the SFR from bursts. 

In the top right panel of Fig. [T?] we show the com- 
plementary quantity pstar, the integrated cosmic stellar 
mass density. All of the C-ACDM models predict a sig- 
nificantly earlier assembly of stars in galaxies than ob- 
servations of high redshift galaxies indicate. The mass 
density of long-lived stars in both our fiducial and the 
HQ model is a factor of ~ 3 higher at 2 ~ 2 and a fac- 
tor of ~ 2 higher at 2 ~ 1 than the observations. How- 
ever, the WMAP3 model produces excellent agreement 
with the stellar mass density as a function of redshift. 
We note here that this tension in the model results (i.e. 
that the C-ACDM model fits the SFR history data bet- 
ter, while the WMAP3 model provides a better fit to the 
stellar mass density) is connected with a possible incon- 
sistency between the two data sets that has b een noted 
recently in several papers (e.g. iHopkins fc Beacom 20061: 
iFardal et al.ll2007l : (Wilkins et al.ll2008l : iDavelbOol ). One 
possible resolution of this tension can be obtained if the 
stellar IMF has changed with time, and was more top- 
heavy at high redshift (we return to this issue in the 
Discussion). 



3.5.3 Evolution of Cold Gas and Metals 

In the bottom left panel of Fig. [14] we investigate an- 
other complementary quantity, the mass density of gas 
that has cooled but not yet formed stars. In our mod- 
els, all of this cold gas is assumed to reside in galac- 
tic disks. One can compare the model predictions with 
the total mass density of Hi gas from blind Hi surveys 
at 2 ~ l|Zwaan et ahlfioosl l. as well as with estimates 
of cold ga s at hi gh redshift from Damped Lyman-a sys- 
tems (e.g. [Prochaska ct al. 2005, and references therein). 
Note that all the observational estimates shown here are 
for atomic gas only, and do not include the contribu- 
tion from molecular gas. Therefore, these observations 
are lower limits for the model predictions of cold gas. The 
C-ACDM models are consistent with the observations to 
2 ~ 4 , but are somewhat low at 2 ~ 5 (note however 



llKomatsu et al.ll2008l) . The new estimates of the parameters 
most relevant to our results, eg = 0.817 it 0.026 and Us = 
0.960^q'Jj3, are intermediate between the values adopted in 
our C-ACDM and WMAP3 models, though somewhat closer 
to the WMAP3 values. 



that iProchaska et al.l l|2005h "caution the reader that the 
results at 2 > 4 should be confirmed by higher resolution 
observations"). However, f2coid in the WMAP3 model is 
much lower than the observations at 2^3, by about a 
factor of two at 2 = 3.5 and an order of magnitude at 
2 ~ 5. 

In the last (bottom right) panel of Fig. Owe show 
the stellar-mass weighted, globally averaged metallicity 
as a function of redshift predicted by the no AGN FB 
model, the HQ model, and the isothermal Bondi models. 
We can compare this with the resu lts of the o bservational 
analysis of SDSS galaxies bV iCallazzi et al. 

I (2007. ) ■ who 
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find a mass-weighted average stellar metallicity at z ~ 
of about solar. The no AGN FB model overshoots this 
value, while the models w i th AG N FB slightly under- 
estimate it. iGallazzi et J] (|2007D showed that the Mil- 
le nnium simulat i ons, b ased on the semi-analytic model 
of ICroton et"ai] l|2006h . give very similar results to our 
models — they underproduce the stellar metallicity at 
2 = by about 20-30 percent. The models predict early 
enrichment, with the mean stellar mass at z 6 about 
25% of solar, and 50% of solar at z ~ 2.5. Because of the 
difficulty of obtaining an unbiased stellar mass weighted 
global mean metallicity at high redshift, we do not at- 
tempt a quantitative comparison with observations, but 
qualitatively, these results seem consistent with the rela- 
tivel y high m etallicities detected in high redshift galaxies 
fe.g. lErb et~ al. 2006). As expected, the WMAP3 model 
predicts somewhat lower metallicity at high redshift. 



3.5.4 The Baryon Budget and its Evolution 

In Fig. [TJl we provide the answer to the question "where 
are the baryons" in our fiducial isothermal Bondi model, 
in both the C-ACDM and WMAP3 cosmologies. AU 
quantities in Fig.[T5]are shown in units of the critical den- 
sity. We show the mass density in collapsed dark matter 
halos above our mass resolution of Mvir ~ 10^" M© mul- 
tiplied by the universal baryon fraction, which represents 
all of the baryons that are available for cooling at a given 
redshift in our model. We also show the mass density of 
hot gas in dark matter halos, and of warm/hot diffuse gas 
that has either been prevented from collapsing into halos 
by the photoionizing background or has been ejected by 
supernova feedback. These two components (warm/hot 
diffuse gas and hot gas in halos) dominate the baryon 
budget at all redshifts, in agr eement with obs ervational 
constraints at low redshift (e.g. lFukugita fc Pe ebles 2004) 
and t he predictions of num erical hydrodynamic simula- 
tions l|Bertone et al.l [20081 . and references therein). Fi- 
nally, we show the baryons in cold gas in galactic disks, in 
stars in the disks and bulges of galaxies, and in "diffuse 
stellar halos" (DSH) around galaxies. Cold gas dominates 
over stars until z ~ 2, where stars begin to dominate. Our 
model predicts that stars in DSH comprise about fifteen 
percent of the mass of stars in galaxies at z ~ 0. 

From this plot we can see that, in the WMAP3 
model, there is just enough baryonic material in collapsed 
halos that can cool at z ;i, 3 to account for the DLAS 
data. If some of the gas that is ejected from galaxies in 
our model instead remains in galactic disks in the form 
of cold gas, perhaps this model could be reconciled with 
the DLAS data. Note, however, that there is no room for 
more cold gas in galaxies at z = 0, so in order to solve 
the problem, the gas needs to be retained at high red- 
shift but ejected or prevented from cooling at low red- 
shift. One can think of plausible physical reasons that 
supernova-driven winds might have more difficulty escap- 
ing galaxies at high redshift, for example, the winds might 
stall against the higher-density IGM that is presumably 
present at these early epochs. 



4 DISCUSSION AND CONCLUSIONS 

We have presented a new semi-analytic model for the self- 
consistent evolution of galaxies, black holes and AGN in 
the framework of the Cold Dark Matter model of struc- 
ture formation. Our models are built on those described 
in SP99, SPF2001, and subsequent papers, but we present 
several important new ingredients here: 

(i) Improved modelling of tidal stripping and destruction 
of orbiting sub-halos and of dynamical friction and satel- 
lite merger timescales. 

(ii) Improved modelling of disk sizes, including realistic 
DM halo profiles and the effect of "adiabatic contrac- 
tion" . 

(iii) A more realistic recipe for star formation in quies- 
cent disks, based on the empirical Kennicutt Law and 
including a surface density threshold for star formation. 

(iv) Updated modelling of the efficiency and timescale 
of starbursts, based on an extensive suite of numerical 
hydrodynamic simulations of galaxy mergers. 

(v) Tracking of a "diffuse stellar halo" component, which 
is built up of tidally destroyed satellites and stars scat- 
tered in mergers. 

(vi) A self-regulated model for black hole growth and 
"bright mode" accretion "lightcurves" based on numeri- 
cal merger simulations with AGN feedback. 

(vii) Galaxy-scale AGN-driven winds, based on numerical 
merger simulations. 

(viii) Fueling of black holes with hot gas via Bondi accre- 
tion, regulated according to the isothermal cooling flow 
model of Nulscn & Fabian ( 200(J). 

(ix) Heating by radio jets, calibrated against observations 
of X-ray cavities in cooling flow clusters. 

We explored the predictions of our new models for 
a broad range of physical properties of galaxies, groups, 
and clusters at z; ~ 0. One key result is an explanation of 
the origin of the characteristic shape of the galaxy stel- 
lar mass or luminosity function. We cast this in terms 
of the "star formation efficiency function" (i.e., the frac- 
tion of baryons in stars as a function of host halo mass) , 
which can be derived empirically by mapping between 
(sub)-halo mass and stellar mass suc h that the observed 
stellar mass funct ion is reproduced (|Moster et al.|[2008l : 
IWang et al.l I2OO6I ). We found that, in our models, the 
shape of this function arises from the fact that supernova- 
driven winds can more efficiently heat and expell gas in 
lower mass halos, while radio-mode heating by AGN is 
more efficient in higher mass halos, because these halos 
have larger mass black holes and these black holes can ac- 
crete hot gas more efficiently. The peak in this function 
at Afvir 10^^ Me occurs because these halos are too 
massive for supernova- driven winds to escape easily, and 
do not have massive enough black holes or a large enough 
virial temperature to fuel efficient radio-mode heating. 

We showed that our model also reproduces the stel- 
lar mass function as a function of galaxy morphology 
(at least for galaxies more massive than lO^^M©), 
the cold gas fractions of disk galaxies as a function of 
stellar mass and the cold gas mass function, the stellar 
mass-metallicity relation for galaxies, and the BH mass 
vs. spheroid mass relation. 

Our model for heating by radio jets is similar in con- 
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Figure 15. The global density in units of the critial density of various components as a function of redshift predicted by the 
isothermal Bondi model (Left: Concordance ACDM with /scatter = 0.4; Right: WMAP3 with /scatter = 0.2). From the highest 
to the lowest curve at 2 = 0, we show: the universal baryon fraction times the mass contained in virialized halos above our 
mass resolution (lO^^M©; solid green); shock-heated hot gas in halos (dotted red); gas in the "intergalactic medium" (see text; 
triple-dot-dashed magenta); stars in galaxies (short dashed orange); cold gas in galaxies (dot-dashed blue); and stars in Diffuse 
Stellar Halos (DSH) around galaxies (long dashed orange). We repeat the observational estimates of Hcoid from Fig. 1141 (squares 
and crosses). 



cept, but different in detailed implementation, to other 
models that have been presented in the literature. We 
tested the basic assumption of our model for fueling of 
BH by h ot gas, the isothe r mal c oolinR flow model pro- 
posed by iNulsen fc FabianI (|2000l ). using observations of 
central temperatures and densities in hot X-ray emit- 
ting gas in nine systems bv lAUen et aLll|200d ). and found 
consistency. We further compared the jet power required 
to solve the overcooling/massive galaxy problem in our 
models with direct measurements of jet power from the 
energetics of bubbles detected i n X-ray gas around ellip- 
tical galaxies (|Allen et al ] |2006l : lRaffertv et aLll2006l ). We 
found tha t our required jet powers lie above the obser- 
vations of lAUen et al. (|2006|) but agree with the higher 
values obtained bv lRafferty et al] l|2006t ). 

We compared the results of our fiducial model with a 
very simple implementation of AGN heating, in which we 
simply switched off cooling in halos more massive than 
~ W^^Mq (the "Halo Quenching" or HQ model). We 
found that this model produced very similar results to 
our fiducial model for the z = stellar mass fraction as 
a function of halo mass (/star(Mhaio)), the stellar mass 
function, the cold gas fractions and cold gas mass func- 
tion, and the global star formation and stellar mass as- 
sembly histories. 

We investigated whether the same models which pro- 
vided a good match to the stellar mass function also 
reproduced the distribution of (specific) star formation 
rates as a function of stellar mass. We found that our 
fiducial model qualitatively reproduces the main features 
of the observed distribution of SSFR vs. stellar mass: 
the models produce a star forming sequence, which dom- 
inates at low stellar masses ( ^ 2 — 3 x IO^^Mq) and a 



quenched population at high stellar masses. The transi- 
tion mass is naturally predicted by the model, and again 
corresponds to the halo mass scale where AGN heating 
becomes effective. This is in marked contrast to models 
without AGN feedback, in which all massive galaxies are 
actively star-forming. 

However, a more detailed comparison with the obser- 
vations indicates that the low-mass star-forming sequence 
in the models occurs at too low a specific star formation 
rate, and does not have the right slope (our star forming 
sequence is flat in SSFR vs. stellar mass space, rather 
than tilted such that low mass galaxies have high SSFR, 
as the observations indicate). This problem is present in 
all of our models, and is independent of AGN feedback 
and robust to the star formation recipe and the values of 
our free parameters. It may be a symptom of the same 
malady that is responsible for producing low-mass galax- 
ies with older ages than those estimated from absorption 
lines in the spectra of nearby galaxies. 

In the Halo Quenching model, the quenching is 
clearly too sharp a function of stellar mass. Essentially 
all galaxies more massive than mstar ~ 10^ ^Mq are com- 
pletely quenched, while the observations indicate that 
massive galaxies have a wide range of specific star for- 
mation rates, from small to moderate amounts of recent 
star formation to activity levels that place them on the 
star-forming sequence. 

We investigated the cosmic histories of the major 
baryonic components of the Universe as predicted in 
our models: star formation and stellar mass, cold gas, 
warm/hot gas, and metals. We found that our fiducial 
Concordance ACDM model produces very good agree- 
ment with the global star formation rate density at z ^ 2, 
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but predicts a higher and flatter SF history at 2 



6 



than most observations indicate. We found that our pre- 
diction of the global SFR density is apparently not af- 
fected by AGN feedback at redshifts above 2^3, because 
most star formation is taking place in small mass halos. 

All of the C-ACDM models predict a significantly 
larger amount of mass in long-lived stars in galaxies at 
redshifts z ^ 0.5, by about a factor of two at z ~ 1 and 
a factor of three at z ~ 2 for the fiducial and HQ mod- 
els. The fact that we reproduce the stellar mass density 
at 2 ~ (by construction) then implies that there is 
not enough evolution in the galaxy population between 
z ^ 1 and 2 ~ in our models with AGN feedback. 
This problem does not seem to be specific to our models, 
but is common to all the recently published models in 
which AGN feedback is implemented in a similar way 
JCattaneo et all l2006l : ICroton et al] l2006l : iBower et all 
|2006| '1. Indeed, one can see from the predictions of the 
global SFR history presented in these works (e.g. Fig. 5 
of Croton et al. 2006, Fig. 8 of Bower et al. 2006) that 
these models all produce very similar predictions for the 
global SF history and for the integrated stellar mass den- 
sity. 

The WMAP3 model, in which structure formation 
occurs later due to the reduced power on small scales, has 
a much lower and more steeply declining global star for- 
mation rate at 2 ^ 2 (about an order of magnitude lower 
than C-ACDM at 2 = 6). The SFR is lower than the 
observational compilation of iHopk ins fc BeacomI (|2006l ) 
at these redshifts by about 0.2-0.3 dex, but still well 
within the observational errors. In contrast to the C- 
ACDM model, the stellar mass density predicted by the 
WMAP3 model is in excellent agreement with observa- 
tional estimates from 2 ^ 4-0. The delayed SF and stellar 
mass assembly history in WMAP3 relat i ve to C-ACDM 
has also been illustrated bv lWang et al.l (|2007t ). 

We also compared the predicted global mass density 
of cold gas in galactic disks Q, cold a, function of redshift 
in our models with estimates of Hi gas mass density at 
2^2;^ 4.5 from Quasar Absorption Systems (Damped 
Lyman-a systems and Lyman limit systems). We found 
that our C-ACDM models had no difficulty producing 
enough cold gas in disks up to 2 4, where the observa- 
tions are the most secure, but the WMAP3 model did not 
fare so well. The predicted values of Slcoid in this model 
were too low by a factor of 2-10 at 2 ^ 3. 

We analyze the redshift-dependent breakdown of 
all the baryons in our fiducial model into each of the 
various components that we track: hot gas in halos, 
warm/hot diffuse gas in the IGM, cold gas in galactic 
disks, stars in galaxies, and stars in Diffuse Stellar Ha- 
los (DSH). We find that hot gas in halos and warm/hot 
gas in the IGM dominate at all redshifts, in agreement 
with t he predictions of nu merical cosmological simula- 
tions (jBertone et al.l |2008| . and references therein) and 
with the observational b aryon census at low redshift 
l|FukuEita fc Peebles! [20oi ) . Therefore, in order to pro- 
duce more cold gas at high redshift, we either require 
more efficient cooling of hot gas or less efficient reheating 
and ejection of cold gas by supernova feedback. Since our 
models are already overproducing stars at high redshift, 
the latter is probably a more promising solution. Our 



models predict efficient early metal enrichment, with the 
average stellar mass- weighted stellar metallicity reaching 
25% of solar at 2 ~ 6 and 50% of solar a.t z ^ 2.5. 

The picture of the cosmic build-up of stars that 
we see in our models is interesting in the con- 
text of a problem pointed out in several recent pa- 
pers (e.g. [H opkins fc Bcacom' l2006l : iFardal et al] l2007l : 
IWilkins etal. 2008; Dave 2008| ): when the best available 
observational estimate of the dust- and incompleteness 
corrected SFR density is integrated, accounting for gas 
recycling under the assumption of a universal stellar IMF, 
the mass of long-lived stars is overestimated by about a 
factor of 2-3. These authors suggest that a possible reso- 
lution to this problem is a non-universal IMF, which was 
more top-heavy at high redshift. In the context of the two 
models we have considered, the C-ACDM models predict 
early structure formation, accompanied by a lot of early 
star formation. We might be able to reconcile these mod- 
els with all the data if in fact the IMF was top-heavy 
at early times, so that most of the star formation that 
we see does not produce long-lived stars. On the other 
hand, the WMAP3 model implies later structure forma- 
tion, and hence less star formation at high redshift. This 
provides an alternate means of reconciling observations of 
the SF history and stellar mass assembly history, which 
requires only that the current observational estimates of 
the SFR at 2 ^ 2 are too high by about a factor of two 
to three. While the WMAP3 picture seems more attrac- 
tive in many respects, it is a concern that the WMAP3 
models do not seem to be able to account for DLAS at 
z again as a consequence of the reduced small scale 
power. Better constraints on the amount of cold gas at 
high redshift from new facilities such as ALMA could help 
to resolve this question. 

In agreement with previous work, we have shown 
that the inclusion of AGN feedback in semi-analytic mod- 
els can plausibly solve the over-cooling problem, the mas- 
sive galaxy problem, and the star formation quench- 
ing problem in the local Universe — a huge step for- 
ward. However, we have also shown that several poten- 
tially serious discrepancies still remain, and we have ar- 
gued that these discrepancies are not peculiar to our 
implementation, but are common to all of the CDM- 
based semi-analytic models currently on the market. 
These discrepancies are connected with low-mass galax- 
ies (m« ^ few X IO^^Mq), however, in the picture that we 
are developing, small galaxies grow into massive gala^xies, 
and the growth of galaxies and AGN are intimately in- 
terconnected, so these problems on small-scales may in- 
dicate or cause more pervasive problems. It is likely that 
these problems are connected to the modelling of cooling, 
star formation and/or supernova feedback, or possibly to 
the CDM power spectrum on small scales. Direct AGN 
feedbaciO is probably not the solution. 

In this paper, we have focussed on predictions of the 
physical properties of galaxies at 2 ~ and the global 



^ By this we mean feedback by an AGN within the galaxy 
itself. Indirect feedback from AGN in external galaxies, 
e.g. via pre- hcating, may be a promising solution (e.g. 
IScannapieco~fc Oh..2004j ) 
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evolution of the major baryonic components of the Uni- 
verse over time. In a planned series of papers, we will in- 
vestigate the predictions of the models we have presented 
hero for multi- wavelength, observable properties (e.g. UV 
through FIR luminosities and colors) of galaxies at low 
and high redsliift, and examine in more detail the distri- 
bution functions and scaling relations of intrinsic galaxy 
properties (e.g. stellar mass, star formation rate, metal- 
licity, etc) at high rcdshift. In addition, wo will explore 
the predictions of our models for AGN properties as a 
function of redshift and environment, and the relation- 
ship between AGN and their host galaxies. 
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